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ABSTRACT 

The  routing  problem  in  a computer-communication  network  is 
modeled  as  a multicommodity  flow  problem.  Two  generalizations  of 
the  EF  algorithm  of  Cantor  and  Gerla  for  solving  the  nonlinear  multi- 
commodity  flow  problem  are  presented.  These  generalizations  are 
of  interest  because  they  retain  the  identity  of  the  components  of  the 
solution  in  an  implicit  and  compact  manner. 

A mini-max  routing  problem  is  defined  as  a linear  program, 
and  a methodology  for  its  solution  as  both  a linear  and  nonlinear 
program  is  presented.  The  two  approaches  are  compared  with  each 
other,  and  a lower  bound  on  the  solution  is  developed  for  the  nonlinear 
approach. 

Two  shortest  route  algorithms  which  are  very  efficient  for 
classes  of  topologies  likely  to  be  considered  in  routing  problems  are 
introduced.  The  shortest  route  problem  is  a subproblem  of  all  of  the 
routing  algorithms  mentioned  above. 

The  computer  load  sharing  problem  is  formulated  as  a multi- 
commodity  flow  problem.  An  algorithm  is  presented  to  solve  the 
load  sharing  problem  which  exploits  the  special  structure  of  the 
problem. 
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CHAPTER  I 
INTRODUCTION 

The  purpose  of  this  thesis  is  to  investigate  and  compare  some  of 
the  algorithms  which  may  be  used  to  find  routing  strategies  in  a com- 
puter-communication network.  Under  a variety  of  assumptions  first 
proposed  by  Kleinrock  [KLEI64]  which  allow  analytic  modeling  of 
packet  switched  computer-communication  networks,  the  static  routing 
problem  can  be  formulated  as  a nonlinear  multicommodity  flow 
problem.  By  the  mechanisms  used  to  solve  such  problems,  one  can 
find  routing  strategies  which  minimize  an  objective  function  such  as 
expected  delay  through  the  network  a randomly  selected  packet. 

The  form  of  a routing  strategy  is  a set  of  parameters  which  specify 
the  fraction  of  packets  arriv'ng  at  each  node  which  is  to  be  sent 
through  each  outgoing  channel  from  the  node  given  the  packet  destina- 
tion. 

Chapter  II  investigates  the  minimum  delay  problem  in  packet 

4 

switched  computer-communication  networks.  Emphasis  is  on  the  EF 
(Extremal  Flow)  algorithm  of  Cantor  and  Gerla  [CANT74]  and  general- 
izations of  that  algorithm.  The  generalizations  of  the  EF  algorithm 
are  interesting  because  they  convey  the  detail  of  the  solution,  i.e.  the 
actual  makeup  of  the  solution  from  which  a routing  strategy  is 
derived,  in  an  explicit  and  compact  way.  The  chapter  discusses  the 
complexities  and  memory  requirements  of  the  various  algorithms  in  a 
comparative  manner. 
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Chapter  III  introduces  two  shortest  route  algorithms  which  are 
efficient  for  certain  interesting  classes  of  networks.  The  shortest 
route  problem  is  an  important  subproblem  of  many  network  problems, 
e.g.  the  EF  algorithm. 

Chapter  IV  discusses  a "throughput"  problem  in  capacitated 
networks.  The  problem  can  be  solved  exactly  via  a linear  program- 
ming approach;  however,  the  linear  program  can  require  a great 
deal  of  computation  time  on  some  relatively  simple  networks.  An 
alternative  is  to  approximate  the  solution  through  a sequence  of 
solutions  to  nonlinear  multicommodity  flow  problems.  The  trade- 
offs between  these  two  approaches  are  discussed. 

Chapter  V is  concerned  with  the  computer  load  sharing  problem. 
The  load  sharing  problem  arises  when  a set  of  computers  v/hich  are 
serving  local  users  in  a batch  processing  mode  are  interconnected 
via  communication  channels  to  share  their  resources  for  better  over- 
all usage  of  the  facilities.  The  algorithm  of  that  chapter  provides 
the  routing  strategy  which  minimizes  the  expected  time  for  passage 
through  such  a system  for  a randomly  selected  job.  The  algorithm 
is  a special  case  of  those  of  Chapter  II;  however,  various  aspects  of 
the  problem  are  exploited  to  increase  the  efficiency  of  the  algorithm. 

The  remainder  of  this  introductory  chapter  is  primarily 
concerned  with  the  nature  of  multicommodity  flows  and  the  solution  of 
linear  multicommodity  flow  problems  in  the  absence  of  capacity 
constraints.  These  concepts  enhance  the  understanding  of  the  decom- 
position methods  for  solving  nonlinear  multicommodity  flow  problems 
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and  for  solving  linear  multicommodity  flow  problems  with  capacity 
constraints.  The  final  section  of  this  chapter  deals  with  the  format 
of  routing  policies,  since  an  optimal  routing  policy  is  the  goal  of  all 
routing  algorithms. 

1.1  Multicommodity  Flows 

The  multicommodity  flow  problem  will  be  stated  in  this  section 
in  the  terms  which  will  be  most  useful  in  later  chapters  for  presenting 
the  body  of  this  report.  The  networks  over  which  the  problems  are 
defined  are  assumed  to  be  directed  and  connected.  Define: 


r. . = 

ij 

value  of  flow  to  be  sent  from  node 

node  j (one  sets  r„  = 0) 

M."* 

II 

flow  on  arc  i 

N 

number  of  nodes  in  the  network 

NA 

number  of  arcs  in  the  network 

C. 

l 

capacity  of  arc  i 

In  the  above,  arcs  are  denoted  by  a single  index,  but  sometimes  it  will 
be  convenient  to  use  the  notation  arc(i,  j)  for  the  directed  arc  from 
node  i to  node  j.  If  the  origin  and  destination  of  arc  k are  OR(k)  = i 
and  DS(k)  = j,  respectively,  then  arc  k corresponds  to  arc(i,  j).  Arcs 
which  have  no  capacity  are  assumed  not  to  exist,  so  that  NA  represents 
the  number  of  arcs  with  capacity  greater  than  zero. 

There  are  three  basic  types  of  constraints  in  dealing  with  net- 
work flows;  nonnegativity  of  flows,  conservation  of  flows,  and  upper 
bounds  on  flows.  The  remainder  of  this  chapter  deals  with  only  the 


( 


-10- 


first  two  types  of  constraints,  and  it  is  assumed  that  if  the  capacity 
of  an  arc  is  greater  than  zero  then  the  flow  through  the  corresponding 
arc  is  allowed  to  be  unbounded.  In  later  chapters,  the  capacity  con- 
straints will  be  included  implicitly  or  explicitly  as  required. 

Writing  the  conservation  equations  and  nonnegativitv  constraints 
for  a general  directed  network  is  a difficult  problem  unless  one 
employs  commodity  distinction  in  the  constraints  (as  will  be  done  later 

■ 

in  this  section).  In  fact,  a systematic  method  for  writing  the  required 
equations  probably  has  not  been  documented.  Some  of  the  conserva- 
tion constraints  may  be  written  in  this  manner:  (i)  choose  anv  set  of 

r ‘1 

arcs  and  denote  this  set  A,  (ii)  find  the  set  B of  origin-destination 
pairs  that  no  longer  have  a path  between  them  if  the  set  A is  removed 
from  the  network,  and  (iii)  write  the  corresponding  conservation 
equation: 

i 

One  makes  the  assumptions  that  all  conservation  equations  mav  be 
written  (and  that  there  are  P of  them),  that  they  are  linear  constraints, 
and  that  they  may  be  written  along  with  the  nonnegativity  constraints  as 

Mf  a b,  f a 0 (1.1.1) 

where  f = [1^,  • • • , b = f b ^ , b9,  • • • , bp]',  and  M has  dimension 

P x NA.  Vector  inequalities  are  interpreted  in  a fashion  similar  to 
the  following  throughout  this  thesis;  f^O  means  f.^0  for  every 
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i «[l,  2,  • • • , NAj . The  set  of  all  flows  f that  satisfy  1.1.1  is 
referred  to  as  the  set  F. 

The  reader  should  not  be  disturbed  by  this  rather  imprecise 
method  of  definition  of  F.  The  development  of  F in  this  manner  has 
had  these  objectives:  (i)  to  point  out  that  there  is  no  existing  system- 
atic method  of  generating  the  constraints  for  F in  a format  which 
does  not  introduce  variables  other  than  those  representing  total  flows 
on  the  arcs,  in  the  hope  that  someone  will  be  inspired  to  solve  that 
problem,  and  (ii)  to  provide  convenient  notation.  It  is  well  estab- 
lished [ASSA75]  that  the  set  F is  a convex  polyhedron  which  implies 
the  existence  of  a constraint  equation  such  as  1.1.1, 


Any  questions  as  to  the  definition  of  F are  answered  by  the  more 
precise  definitions  of  F provided  by  the  following  two  constraint 
formats.  They  allow  a compact  representation  of  the  constraints; 
however,  commodity  distinction  is  required. 

In  the  N-commodity  formulation,  a commodity  is  the  ensemble 
of  inputs  to  the  network  which  are  destined  to  a particular  node.  The 
conservation  equations  are  written  at  each  node  j for  each  commodity 


k as 


i such  that  i such  that 
OR(i)  = j DS(i)  = j 


where  fv  is  the  flow  in  arc  i which  is  destined  for  node  k.  Collecting 
the  equations  for  each  commodity,  employing  matrix  notation,  and 
adding  the  nonnegativity  constraints,  the  result  is  denoted 


r 
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Ef14  = 


P * 0 


(1.1.2) 


where  E has  dimension  N x \TA  and  its  elements  are  tl  or  0,  = 

^2*  * * * » ’ anc*  ^ = . bg , * • • , b^]',  where  tr  = - L if  j = 

k ^ 

and  b.  = r^  if  j r k.  E = fe. .]  is  generally  called  the  node-arc 

incidence  matrix  and  is  independent  of  the  destination  node,  k,  and 

6ij  = ^ ^ OR(j)  = i»  e.^  = -1  if  DS(j)  = i,  and  e..  = 0 otherwise.  The 

set  F is  also  defined  by  1.  1.2;  in  particular,  fcF  if  and  only  if  there 

exist  f , r,  * • • , f such  that  f = L / and  r satisfies  1.1.2  for  ever-. 

k = 1 

ke{  1, 2,  • • •,  n)  [ASSA75] . 

O 9 

The  N -commodity  formulation  (actually  there  are  % - X 
commodities)  is  closely  related  to  the  N-commoditv  formulation 
except  that  a commodity  is  associated  with  both  an  origin  and  a 
destination.  The  conservation  equations  are  written  at  each  node  j 
for  each  commodity  as 


mk 


if  j = m 


T jjnk 

L i 

i such  that  i such  that 


-r  if  i - k 

mk 

0 otherwise 


I Ck 

such  th; 

OR(i)  = j DS(i)  = j 

where  is  the  flow  in  arc  i with  origin  rp  and  destination  k. 
Collecting  the  equations  for  each  commodity,  employing  matrix  nota- 
tion, and  adding  the  nonnegativity  constraints,  the  result  is  denoted 


Ef 


mk 


mk 


,JTik 
1 2 


(1.1.3) 


where  E is  the  same  matrix  as  in  1.1.2,  fmk  = [f"',k,  fink>  • • • , f™k]' 

1 *- 


umk  rvmk  ...  imk„  u .mk 

and  b = |b j , 02  , * * * , b^.  j' , where  b = r 


mk 


if  j = m,  b‘ 


mk 
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-r  . if  j = k,  and  b.  =0  otherwise, 
mk  J ] 


1.2  Extremal  Flows,  Trees,  and  Chains 


This  section  will  investigate  the  nature  of  solutions  to  linear 

programs  when  the  constraints  have  the  form  of  1.1.1,  1.1,2,  or 

1.1.3.  Depending  on  the  particular  format,  solutions  wrill  be  called 

extremal  flows,  trees,  or  chains.  The  powerful  decomposition 

methods  of  later  chapters  will  be  explained  in  terms  of  these  notions. 

NA 

The  problem  now  under  consideration  is  to  minimize  L c.f.  A «_*'f 

i=l  11  = 

such  that  fcF.  This  linear  program,  also  called  a shortest  route 
problem,  is  a subproblem  of  the  algorithms  which  employ  the  decom- 
position methods.  The  method  of  attack  will  be  to  solve  the  problem 
subject  to  1.1.3  for  each  of  the  commodoties;  superposition  of 
appropriate  solutions  yields  the  N-commodity  solutions,  which  in  turn 
are  superimposed  to  yield  the  solution  without  commodity  distinction. 
Consider  the  following  problem: 

min  c'f  such  that  Efmk  = bmk  , fmk  * 0 (1.2.  1) 


Generally  the  solution  is  not  bounded  since  the  set  F is  not  bounded. 

However,  in  the  case  that  c ^ 0 (this  assumption  is  implicit  in  the 

remainder  of  this  section),  then  a bounded  solution  to  1.2.  1 does  exist. 

Assume  that  commodity  (m,k)  obeys  a pure  routing  strategy  which  can 

be  described  by  a chain,  (m  = nQ,  n^,  • • • , m - k),  or  as  a sequence  of 

arcs,  (i  , i„,  • • • , i .),  such  that  OR(i  ) = n , and  DS(i  ) = n for 

1 2 j P p-1  P P , 

p<£l , 2,  • • • , j} . The  corresponding  flow,  f£nk  = (rmk  ^ ^1*  *2.  * * ’ 

1 0 otherwise 

clearly  satisfies  1.  1.3,  and  the  cost  of  such  a solution  is 
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r ( L Ci  ).  If  c.  is  thought  of  as  the  distance  of  arc  i,  and  the 
mK  p=l  T)  1 

length  of  a chain  is  the  sum  of  the  distances  of  its  corresponding  arcs, 
then  minimization  of  the  cost  of  a pure  routing  strategy  from  node  m 
to  node  k can  be  seen  as  equivalent  to  finding  the  minimum  length 
chain  from  node  m to  node  k.  If  the  minimum  length  chain  is  unique, 
then  the  corresponding  flow,  denoted  O^'  , is  the  unique  solution  to 
1.2.1.  If  there  are  j minimum  length  chains,  whose  corresponding 
flow  vectors  are  denoted  c^1^,  • • then  any  convex  combination. 


f111^,  of  the  chains  is  a solution  to  1.2.  1;  for  this  case,  fm^  = ^ 

i=l  1 l1 


mk 

) 


such  that  X.  ^ 0 and  L = 1. 
l i = 1 i 

The  ambiguhv  that  arises  from  degeneracy,  i.  e.  more  than  one 

minimum  length  path,  needs  to  be  resolved.  Later  in  this  thesis, 

reference  will  be  made  to  the  set  of  chains  which  solve  1.2.  1 for  a 

given  metric  c.  This  set  contains  exactly  one  chain  for  each 

commodity.  Degeneracy  is  allowed  to  be  resolved  arbitrarily  as  long 

as  one  is  consistent.  Consistency  in  this  matter  means  that  if  the 

declared  shortest  chain  between  node  m and  node  k is  (m  = n . n . 

o 1 

• • • , m =k),  then  the  declared  shortest  chain  between  node  q and  node  k 
where  q = n for  pcCl,  2,  • • • , i-lJ  must  be  (q  = n , n n.  = k). 

p ^ J n p p- 1 i 

At  this  point,  one  should  be  aware  of  the  reasons  that  the 
solution  to  1.2.1  is  not  bounded  in  general.  If  the  length  of  a closed 
path  is  negative,  then  an  arbitrarily  large  flow  may  circulate  or  that 
closed  path  while  still  satisfying  1.  1.1.  Therefore,  in  this  case,  the 
objective,  c'f,  has  no  finite  lower  bound.  However,  by  specifying  that 
the  metric,  c,  must  be  nonnegative,  the  existence  of  a bounded  solution 


f 
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to  1.2.1  is  assured. 

If  one  has  the  capability  to  solve  1.2.1,  then  the  solution  to  the 
following  problem  is  at  hand: 


minimize  c'f  such  that  Ef^  = b\  f%  ^ 0 


<1.2. 2) 


In  particular,  if  one  has  solved  1.2.1  for  every  m([l,  2,  • • • , x)  - tk] , 

then  the  solutions  may  be  superimposed  to  solve  1.2.2.  The  solution 

to  1.2.2  is  denoted^,  and  $ is  related  to  the  set  of  chains  whicli  solv« 

1.2.1  by  ft  = ^L  <Pm*<.  The  care  in  resolving  degeneracy  is  rr.ani- 
k 

fested  in  the  resultant  form  of<dS  the  solution  to  1.2.2  is  called  a 


tree  because  ft  has  exactly  N- 1 nonzero  elements  when  r > 0 foi 

mk 


every  mcll,  2,  * • • , NJ  such  that  m f-  k,  and  the  nonzer  o element  corres- 
pond to  a set  of  arcs  which  is  a tree  of  the  network.  Each  tree  could 
be  considered  rooted  at  the  destination  node  with  each  arc  directed 
towards  that  node. 

Finally,  one  may  minimize  c'f  subject  to  1.1.1  by  superimposing 
the  set  of  chains  which  solve  1.2.1  or  the  set  of  trees  which  solve 
1.2.2.  The  solution,  <p,  is  called  an  extremal  flow  or  a shortest 


route  flow,  and  a relationship  for  the  extremal  flow  in  terms  of  trees 


N i-  in  in  rnk 
and  chains  is  <P  = L p = 2 2 ^ . 

k = 1 k = 1 m=^l 
mrk 


N N 


1 . 3 The  Representation  of  Routing  Policies 

One  relatively  conpact  wav  to  state  a routing  policy  is  to  provide 


for  each  node,  m,  and  each  possible  destination,  k,  where  k r m,  a 

NA 

set  of  parameters  Ttmk  such  that  2 = 1,  trT11^  ^ 0,  and  = 0 if 

i i=l  i t l 


rrik 

OR(i)  r m.  IT.  is  the  fraction  of  How  arriving  at  node  m which  is  to 
be  sent  next  along  branch  i to  the  node  DS(i).  The  practical 
importance  of  such  a representation,  in  the  context  of  packet  switched 
networks,  is  that  the  policy  requires  a small  amount  of  information 
about  any  given  packet  (i.  e.  the  destination)  in  order  for  the  packet  to 
be  processed  upon  its  arrival  at  any  given  node. 

The  ultimate  goal  of  each  routing  algorithm  is  to  provide  such 
a routing  policy,  although  sometimes  this  is  not  evident.  Routing 
policies  are'often  implicit  in  routing  algorithms,  and  the  desired 
explicit  format  is  quite  difficult  to  extract  from  the  "solution"  to  a 
given  pioblem  when  this  situation  occurs.  These  thoughts  will  be 
given  fuller  development  in  the  chapter  which  follows. 


CHAPTER  II 

NONLINEAR  MULTICOMMODITY  FLOW  PROBLEMS 

This  chapter  deals  with  the  application  of  nonlinear  multi- 
commodity  flow  models  to  the  determination  of  a routing  police  for  a 
computer-communication  network.  Kleinrock's  model  of  a packet 
switched  network  as  a network  of  queues  shall  be  taken  as  the  founda- 
tion of  such  an  application  [KLEI64]. 

The  Model.  In  this  model,  there  exists  a computer  at  each 
node  of  the  network.  The  computers  may  perform  a variety  of 
services  (in  terms  of  the  data  processing  of  packets)  which  are  aimed 
towards  the  errorless  delivery  of  messages  with  minimum  delay 
between  the  network's  subscribers.  A long  message  may  be  broken 
into  smaller  units  which  are  denoted  packets.  The  packets  propagate 
through  the  network  as  separate  entities  and  are  reassembled  at  the 
destination  node  into  the  original  message.  Packets  travel  through 
the  network  in  a store-and-forward  fashion.  That  is,  each  nodal 
computer  observes  the  destination  of  every  arriving  packet.  If  the 
packet  is  destined  elsewhere,  the  computer  selects  an  appropriate  out- 
going channel  and  places  the  packet  in  the  queue  associated  with  that 
channel  to  be  forwarded  in  its  turn. 

The  problems  associated  with  packet  switched  communication 
networks  are  numerous,  and  only  a few  of  these  will  be  mentioned. 
First  is  the  location  of  the  nodes,  the  selection  of  the  topology,  the 
choice  of  the  transmission  medium,  and  the  selection  of  capacities. 
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( Error  control  is  not  only  important  on  each  individual  link,  but  also 

from  source  to  destination,  since  packets  may  be  lost  for  various 
reasons.  Finite  storage  facilities  may  require  reservation  of  memory 
for  the  reassembly  of  messages.  Descriptions  of  the  state  of  the  net- 
work may  be  desired  to  allow  the  network  to  adapt  to  failure  or  heavv 
loading  of  nodes  and/or  links,  and  then  one  must  decide  on  how  and 
when  to  send  the  network  description. 

This  chapter  addresses  only  one  of  these  problems,  namely  the 

routing  problem,  i.e.  how  should  a nodal  computer  select  the  output 

channel  for  any  given  packet  which  may  arrive  at  its  node.  This 

problem  is  especially  important,  however,  for  several  reasons.  It 

is  related  to  the  topology  and  capacity  selection  problem  because,  to 

evaluate  the  relative  merits  of  different  selections,  one  must  specify 
I 

how  the  network  will  be  used  before  one  can  determine  what  demands 
will  be  placed  on  its  various  components.  Moreover,  the  routing 
problem  is  important  for  more  than  just  the  initial  design  of  the  net- 
work. It  is  also  important  in  maintaining  the  effectiveness  of  the  net- 
work during  normal  operating  conditions  and  during  emergency 
conditions  such  as  failure  of  components  of  the  network,  and  in  deter- 
mining how  the  network  should  grow  with  increasing  usage.  To  lose 
traffic  because  the  network  cannot  adapt  to  unusual  conditions  or  to 
invest  capital  into  a larger  network  when  a more  robust  routing 
strategy  would  suffice  are  design  errors  which  unnecessarily  degrade 


the  cost-effectiveness  of  the  network. 

The  simplest  model  of  a packet  switched  network  having  some 
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degree  of  realism  is  a network  <>!'  queues  model.  This  approach  was 
first  developed  by  Kleinrock  [KLUI64],  All  delays  are  assumed  to 
be  negligible  except  queueing  delays  which  arise  from  the  finite  trans- 
mission rates  over  channels  (determined  by  the  channel  capacities) 
and  the  randomness  of  arrivals  of  packets  which  are  to  be  sent  over 
the  various  channels.  One  also  assumes  in  this  model  that  each 
packet  must  have  fully  arrived  at  a node  prior  to  its  forwarding,  and 
that  the  queues  associated  with  the  various  channels  have  infinite 
storage  capability.  Then  the  following  additional  assumptions  or 
approximations  are  made  in  order  to  apply  the  results  of  Burke 
[BURK56]  and  Jackson  [JACK63]  concerning  networks  of  queues, 
making  the  problem  analytically  manageable: 

(i)  The  arrival  times  of  packets  with  an  origin  i and 

destination  j form  a Poisson  process  with  mean  arrival 
rate  r„  (bits/second)  which  is  independent  of  all  other 
network  processes. 

(ii)  The  length  of  a packet,  which  is  proportional  to  the 

service  time  at  a node,  is  an  exponentially  distributed 
random  variable  with  mean  1 In  (bits). 

(iii)  The  length  of  a packet  is  assumed  to  be  chosen 

independently  of  all  network  processes  upon  arrival  at 
each  successive  node  from  source  to  destination. 

Under  all  of  these  assumptions,  the  average  delay  of  a randomly 
selected  packet,  T (seconds /packet),  becomes: 


J 


i 


-20- 


NA  f 

1 Y 1 

Jly  L TT^T 

i=l  1 1 


where  f.  is  the  average  bit  rate  (bits /second)  on  the  ith  channel  and  is 

determined  by  the  demands  and  the  routing  policy,  and  where  y - 

L r . .. 
i,j  bl 

The  expression  for  T can  become  more  elaborate  if  causes  for 
delay  other  than  queueing  are  included  in  the  model,  e.g.  nodal 
processing  time,  propagation  time,  network  control  traffic,  etc. 
[K.LKI70],  The  methods  of  this  chapter  will  remain  valid  as  long  as 
the  delay,  expressed  as  a function  of  How,  T(f),  is  convex  and  the 
gradient  of  T,  VT  = |CT /Of  , 0T /Sf9  , • • ■ , ST/bf^.  ^ ]' , is  nonnegative  over 
the  set  FOG,  where  f € G if  0 £ f.  < C.  for  every  i*[l,  2,  • • • , NAj . 

Previous  Work.  One  criterion  which  is  used  to  determine  an 
'optimal'  routing  policy  is  the  minimization  of  the  expected  delay 
through  the  network.  Such  a routing  policy  is  termed  optimal  in  the 
static  sense  with  respect  to  average  delay.  One  can  update  the  policy 
as  one  acquires  new  estimates  of  the  demands,  i.e.  the  set  of  require- 
ments, {r_i.  In  this  case,  the  algorithms  based  on  the  above  criterion 
may  be  termed  quasistatic,  but  are  not  necessarily  optimal  in  any 
sense,  although  perhaps  near  optimal  in  certain  senses  under  certain 
circumstances. 

The  algorithms  of  this  chapter  are  termed  centralized,  in  the 
sense  that  the  required  computations  are  most  efficiently  performed  at 
a central  location  in  the  network.  This  of  course  entails,  for  an 
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• actual  quasistatic  implementation,  the  transfer  of  some  "state  of  the 

network"  description  to  a central  location,  a calculation  of  a desirable 
routing  policy,  and  then  distribution  of  this  routing  information 
throughout  the  network.  It  seems  possible  that  a distributed  algorithm 
may  reduce  this  communication  problem  at  a small  cost  in  nodal 
computer  complexity.  The  reader  is  referred  to  the  work  of  Gallager 
[GALL77]  for  an  investigation  of  this  possibility,  and  for  other 
references  concerning  distributed  routing  algorithms.  The  choice  of 
distributed  versus  centralized  routing  in  an  actual  implementation 
probably  should  be  based  on  many  factors  such  as  the  reliability  of  the 
components  of  the  network  and  the  nature  of  the  demands  on  the  system. 
Although  information  must  be  passed  farther,  in  general,  for  a 
centralized  procedure,  qualitatively  one  can  see  that  there  may  exist 

( 

cases  where  the  information  can  be  passed  less  often.  A distributed 
algorithm  converges  to  the  solution  of  a particular  static  problem  by 
the  process  of  information  interchange  among  the  nodes.  On  the  other 
hand,  in  the  case  of  a centralized  algorithm,  information  is  passed 
only  for  gathering  data  and  communicating  the  routing  policy  to  the 
nodes. 

« 

The  algorithms  to  which  this  chapter  will  be  directed  are  also 
referred  to  as  exact  rather  than  heuristic.  Exact  in  this  case  means 
that  one  can  come  arbitrarily  close  to  the  solution  of  the  minimum 
delay  problem  using  such  algorithms,  given  an  exact  computing 
machine  and  a sufficient  number  (which  is  finite)  of  iterations.  Exact 


algorithms  appear  to  have  the  same  order  of  complexity  as  heuristic 
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algorithms  exemplified  by  those  of  Frank  and  Chou  [FRAN71a|  and 
Fultz  [FULT72],  Moreover,  exact  algorithms  have  the  advantage  of 
always  having  the  capability  of  being  driven  to  a near  optimal  solution, 
as  opposed  to  heuristic  algorithms  which  invariably  seem  to  perform 
poorly  in  some  cases. 

Furthermore,  this  chapter  is  concerned  with  algorithms  which 
retain  the  identity  of  components  of  the  total  flows  on  each  link  in  the 
solution.  Otherwise,  a routing  policy  is  not  easily  derived.  An 
algorithm  which  performs  the  minimization  without  retaining  this 
identity  is  the  Frank-Wolfe  (FW)  algorithm  [FRANTlb].  A very 
similar  algorithm  was  denoted  the  flow  deviation  (FD)  algorithm  and 
was  initially  applied  to  packet  switched  communication  networks 
[FRAT73]  . The  applications  of  the  FW  algorithm  to  transportation 
networks  were  studied  by  Leblanc,  Morlock,  and  Pierskalla  [LKHL75] 
and  by  Nguyen  [NGUY741.  A variation  on  the  FW  algorithm  wa.-; 
proposed  by  Klessig  [KLES74]  who  altered  the  approach  to  the  one- 
dimensional minimization  of  the  algorithm.  Accordingly,  the  FW 
algorithm  arrives  at  a value  of  minimum  delay  but  does  not  reveal  how 
to  achieve  this  value  through  routing.  One  consequence  of  this  is  the 
inability  to  arbitrarily  perturb  the  demand  matrix,  and  then  arrive  at 
the  new  solution  without  completely  resolving  the  problem  (see 
Appendix  B).  Naturally  the  memory  associated  with  such  algorithms 
is  comparatively  small,  and  the  FW  algorithm  can  handle  quite  large 
networks  more  easily  than  algorithms  that  are  capable  of  detem  ining 
a routing  policy.  The  use  of  the  FW  algorithm  is  directed  towards  a 


criterion  for  comparing  networks  when  the  network  is  in  its  first 
design  stage.  The  FW  algorithm  (or  equivalently  the  FD  algorithm) 
was  encoded  as  a part  of  this  thesis  in  order  to  serve  as  a standard  of 
comparison  for  the  computations  cited  in  Chapter  IV. 

Some  multicommodity  flow  problems  have  methods  which  are 
effective  for  a small  number  of  commodities.  In  particular.  Golden 
developed  an  algorithm  for  import-export  routing  involving  two 
commodities  [GOLD75],  and  Schwartz  and  Cheung  developed  an 
algorithm  suited  to  computer-communication  networks  with  a very 
small  number  of  origin-destination  pairs  [SCHW76],  The  algorithms 
under  consideration  in  this  chapter  are  directed  to  those  cases  where 
the  demand  matrix,  R = [r- ^ ] , has  many  nonzero  terms,  although  they 
should  not  necessarily  perform  poorly  when  that  is  not  the  case. 

The  first  section  of  this  chapter  outlines  the  algorithm  of  Cantor 
and  Gerla  [CANT74]  which  is  called  the  EF  algorithm,  and  the 
remainder  is  devoted  to  related  algorithms  based  on  alternative 
decomposition  schemes.  The  object  of  these  algorithms  is  to  mini- 
mize T(f)  subject  to  feF,  and  to  find  the  associated  routing  policy. 
Notice  that  T(f)  only  makes  sense  when  f satisfies  the  capacity 
constraints,  so  that  T(f)  ^ =>  if  f.  2 C.  for  any  i e Cl , 2,  • • • , NAj . 


2.  1 The  EF  Algorithm 

The  short  presentation  of  this  section  views  the  EF  algorithm  as 
a column  generating  programming  technique  [ASSA75j.  The  name  EF 
is  drawn  from  the  term  "extremal  flow"  which  was  introduced  and 
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defined  in  Section  1.2.  Initially,  one  assumes  that,  a restricted  set  of 

extremal  flows,  called  a basis  ^ j ), ),  • • • , ^ - )J , has  been  provided. 

Let  f be  a convex  combination  of  the  elements  of  the  basis,  i.e. 

f = i such  that  Z X,  = 1 and  * 0.  Clearlv,  fcF;  however, 

i _ 1 i u;  i=  1 1 1 

not  every  element  of  F can  be  expressed  as  a convex  combination  of 
the  given  basis  elements.  The  restricted  master  problem  is 

1 

minimize  T(  Z X-'Ptx) 
i = 1 ! 

such  that  t X = 1,  X-  ^ 0 
i = 1 1 1 

which  can  be  solved  as  a multidimensional  minimization  over  the 

parameters  under  constraint.  Denote  the  solution  to  the 

1 ^ J 

restricted  master  as  f’ , which  is  said  to  be  the  current  estimate  of 
the  solution  of 

minimize  T(f)  such  that  fcF  (2.1.1) 

a. 

Since  F is  a convex  set  and  T(f)  is  a convex  function,  then  fr  is 
the  solution  to  2.1.1  unless  there  exists  fcF  such  that 

Rel  Cost  (f)  £ <VT(f*),  <f-f*)>  < 0 (2.1.2) 

where  (,  ) denotes  the  usual  inner  product  in  Euclidean  space.  This 
condition  for  optimality  is  derived  from  the  Kuhn-Tucker  theorem 
[VARA72]  by  Cantor  and  Gerla  [CANT74]  . Define  c = 7T(f);  then  the 
minimization  of  Rel  Cost  (f)  such  that  fcF  is  equivalent  to  the  minimiza- 
tion of  c'f  subject  to  1.  1.  1.  A solution  to  this  problem  is  an  extremal 
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flow,  which  is  denoted  ^ ^ y and  the  method  of  solution  is  by  solving 

a shortest  route  problem  as  outlined  in  Section  1.2.  This  shortest 


route  calculation  is  called  the  subproblem  of  the  EF  algorithm.  If 


the  solution  to  the  subproblem,  is  such  that  Rel  Cost  . ()  ^ 0, 

then  f solves  2.  1.  1 and  the  algorithm  terminates;  otherwise,  O. 
is  added  to  the  restricted  set  of  extremal  flows,  and  the  master 


problem  is  resolved  with  the  assurance  that  the  objective  function, 
T(f),  can  be  lowered,  and  so  forth.  The  number  of  extremal  flows  is 


finite,  and  if  one  can  solve  the  restricted  master  problem  exactly, 
then  it  follows  that  the  algorithm  solves  the  problem  in  a finite  number 
of  iterations. 


Implicitly,  in  the  above,  the  first  estimate  of  the  solution  to 
2.1.1  is  assumed  to  satisfy  the  capacity  constraints  since  the  multi- 
dimensional minimization  of  the  master  problem  requires  a well 
defined  gradient.  The  starting  problem  is  that  of  finding  such  an 
initial  estimate  of  the  solution;  this  problem  is  discussed  further  in 
Chapter  IV  where  an  easily  implemented  procedure  for  starting  the 
algorithm  is  presented. 

The  analogy  of  the  EF  algorithm  with  linear  programming  and 
the  simplex  method  are  apparent.  The  columns  of  the  program  are  all 
possible  extremal  flows.  Only  a small  number  of  these  are  known 
explicitly,  namely  the  members  of  the  basis,  which  in  turn  are  the  only 
columns  with  which  a nonzero  primal  variable,  can  be  associated. 

At  the  conclusion  of  the  master  problem  (assuming  perfect  minimization), 
the  cost  with  respect  to  the  current  basis,  Rel  Cost  (<P^.  ^),  is  zero  for 
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every  basis  element  such  th.  ■ . > 0 (the  remainder  of  the  basis 
elements  may  be  deleted  in  • EF  algorithm  for  more  efficient  use  of 
memory  in  the  computer).  Corresponding  to  the  simplex  algorithm's 
search  for  the  nonbasic  column  with  minimum  relative  cost,  the  sub- 
problem generates  such  a column,  <p(  ^ In  the  nonlinear  problem, 

Rel  Cost  (<P,.  is  not  only  related  to  the  ability  of  <p, . , . to  lower  the 
(j+1)  “ (j+1) 

objective  function  upon  being  included  in  the  basis,  but  also 

* M * 

- Rel  Cost  (<p^+  is  an  upper  bound  to  T(f  ) - T(f  ),  where  f is  the 
current  estimate  of  the  optimal  solution  found  by  the  master  problem, 

M 

and  f is  the  solution  to  2.  1.  1 (see  Appendix  A,  Thm.  2.  1A). 
Theoretically,  one  stops  the  algorithm  when  Rel  Cost  (<P^  0; 

however,  as  a practical  matter,  the  upper  bound  is  a powerful  tool  in 
knowing  when  stopping  the  algorithm  is  appropriate. 

The  multidimensional  minimization  which  is  required  in  the 
master  problem  was  solved  via  the  gradient  projection  method  by 
Cantor  and  Gerla  [CANT74],  This  technique  reduces  the  problem  to 
a series  of  one-dimensional  minimizations  each  of  which  is  along  a 
feasible  downhill  direction  (but  not  the  direction  of  steepest  descent). 
Best  employed  a variation  of  this  approach  towards  the  multidimensional 
minimization  problem  which  may  be  of  some  benefit  [BEST75] ; however, 
there  is  not  an  actual  comparison  showing  net  reduction  of  overall  com- 
putation time  due  to  his  alternate  approach,  even  though  his  algorithm 
was  actually  encoded. 

This  section  is  concluded  with  some  remarks  on  finding  the 
optimal  routing  strategy  from  the  solution  to  the  EF  algorithm. 


t 
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Knowing  the  optimal  flow  f to  the  problem  2.1.1  does  not  imply  the 
solution  to  the  routing  policy  because  f is  a convex  combination  of 
flows  without  commodity  distinction.  Specification  of  an  optimal 
routing  polic>,  however,  requires  at  least  a knowledge  of  the  destina- 
tions of  tne  flows  making  up  the  solution.  For  example,  if^. ) is  the 

component  of  ' destined  for  node  k,  then  f*  = ^ \.<p  = 1 L v«>k 

w . i k i=1  1 M 

and  r **  .2^  is  the  component  of  the  optimal  flow  destined  for 

node  k.  Consequently,  the  routing  variables  defined  in  Section  1.  3 

may  be  calculated  as  follows; 

fk# 

(2,i,3) 

j 

m such  that 
OR(j)  = m 

Clearly,  the  optimal  routing  policy  is  known  only  if  sufficient  informa- 
tion is  stored  with  each  extremal  flow  to  be  able  to  reduce  it  into  its 
corresponding  trees.  The  most  compact  way  to  retain  this  informa- 
tion (and  the  way  suggested  by  Cantor  and  Gerla)  is  via  a routing 
matrix,  NX,  of  dimension  N x N where  NX„  is  the  next  arc  from  node 
i on  the  shortest  route  to  node  j.  Thus,  having  stored  a routing 

i 

matrix  for  each  extremal  flow,  when  the  EF  algorithm  terminates,  the 
necessary  trees  may  be  constructed  which  allow  solution  of  2.  1.  3 for 
every  pair  of  m and  k. 

2. 2 The  TR  Algorithm 

The  TR  algorithm  is  an  original  variation  on  the  Cantor  and 
Gerla  routing  algorithm  with  a decomposition  based  on  trees  rather 


than  extremal  flows.  The  computational  experience  with  the  TR 
algorithm  which  will  be  presented  in  Chapter  IV  and  the  theoretical 
comparisons  in  the  remainder  of  this  chapter  suggest  that  the  TR 
routing  algorithm  should  be  competitive  in  speed  with  the  EF 
algorithm  in  solving  2.  1.  1,  and  makes  more  efficient  use  of  memory 
in  a computer  implementation.  Moreover,  an  optimal  routing  policy 
is  determined  at  the  conclusion  of  the  algorithm  with  far  less  effort. 

This  is  the  master  problem  for  the  TR  algorithm; 


minimize  T(f)  = TIE  L X,  cp  . ) 

k = l i=l  '-k  (i)/ 


N lk 


% 


such  that  L 

i = 1 lk 


■ '•  * 0 


(2.2.1) 


where  each^.j  is  a member  of  a given  restricted  set  of  trees  (as 
defined  in  Section  1.2)  which  is  called  a basis,  and  t denotes  the 

K * 

number  of  trees  in  the  basis  associated  with  flows  destined  for  node  k. 
In  2.2.  1,  f is  said  to  be  a convex  combination  of  the  trees  in  the  given 
basis.  Denote  the  solution  to  2.2.  1 as  f . 

The  subproblem  to  the  TR  algorithm  is  analogous  to  that  of  the 
ET  algorithm;  find  the  set  of  trees  which  minimizes  c'f  subject  to 

# s{c 

1.1.2  where  c = VT(f  ).  If  f is  not  the  global  optimum,  then  newly 
generated  trees  are  added  to  the  existing  basis  and  the  algorithm  is 
iterated. 

At  first  glance,  this  master  problem  may  look  more  difficult  to 
solve  than  the  master  problem  of  the  EF  algorithm  due  to  the  increase 
of  the  number  of  columns  and  associated  variables.  On  the  other  hand, 
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an  advantage  of  the  format  of  2.2.  1 is  apparent.  Suppose  a restricted 
set  of  trees  was  generated  by  the  same  shortest  route  problems  which 
generated  a restricted  set  of  extremal  flows.  In  this  case,  the  region 
of  F which  can  be  expressed  as  a convex  combination  of  the  trees 
includes  the  region  which  can  be  expressed  as  a convex  combination  of 
the  extremal  flows.  If  the  respective  master  problems  were  solved 
at  this  point,  then  the  TR  master  problem  would  generally  find  a better 
estimate  of  the  solution  to  2.1.1. 

Section  2.3  will  deal  more  precisely  with  the  comparisons  of  the 
master  problems;  the  remainder  of  this  section  compares  memory 
requirements  and  ease  in  reducing  solutions  to  optimal  routing 
policies.  Cantor  and  Gerla  showed  that  the  solution  to  2.  1.  1 can  be 
expressed  as  a convex  combination  of  at  most  NA+1  extremal  flows 
[CANT74].  Via  virtually  identical  arguments,  the  solution  to  2.  1.  1 
can  be  expressed  as  a convex  combination  of  at  most  NA+N  trees  (see 
Appendix  A,  Thm.  2. 2.  A).  Consider  the  storage  implications  of 
these  facts  (it  will  be  assumed  that  both  integers  and  reals  require  one 
memory  location,  although  on  some  machines  integers  may  be  stored 
more  compactly  than  reals).  Trees  have  no  more  than  N-l  nonzero 
terms,  and  if  they  are  stored  along  with  an  appropriate  arc  index  for 
each  term,  then  the  maximum  storage  requirement  is  2(N-1)(NA+X) 
locations  for  the  basis.  The  corresponding  maximum  requirement  for 
the  EF  basis  is  NA(NA4l).  If  the  number  of  arcs  exceeds  three  times 
the  number  of  nodes,  then  the  former  term  is  smaller. 

Now  suppose  that  the  optimal  routing  policy  is  desired  at  the 
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conclusion  of  the  algorithm,  so  that  one  must  also  store  a routing 
matrix  for  each  extremal  point.  The  maximum  memory  requirement 

9 

associated  with  the  EF  basis  is  now  (N“+NA)  (NA+1)  which  is  about  a 
factor  of  N greater  than  the  required  storage  for  the  TR  algorithm. 

In  order  to  find  the  optimal  routing  policy,  one  must  find  f^"  for  every 
k({l , 2 , • • • , nJ.  For  the  TR  algorithm,  this  requires  at  most 
(N-lHNA+N)  multiplications  in  a very  straightforward  procedure. 

After  one  has  found  the  trees  associated  with  all  the  extremal  flows  in 
the  solution  of  the  EF  algorithm,  combining  them  to  form  every 
requires  a maximum  of  (N- 1 )(NA+ 1 )(N)  multiplications  and  additions. 

Evidently,  the  TR  algorithm  looks  very  attractive  in  both 
storage  requirements  and  convenience  of  the  format  of  the  solution  for 
establishing  the  optimal  routing  policy.  Since  the  subproblem  of  the 
EF  algorithm  is  of  exactly  the  same  complexity  as  that  of  the  TR 
algorithm,  only  a comparison  of  the  efficiency  of  the  master  problems 
remains  to  be  examined.  The  next  section  deals  with  a suggested 
method  for  solving  the  multidimensional  minimization  of  the  TR  master 
problem  and  discusses  the  complexity  in  a comparative  manner  with  that 
of  the  EF  master  problem. 

2 . 3 The  TR  Master  Problem 

Several  concepts  presented  in  Section  2.  1 will  be  generalized  in 
order  that  the  TR  algorithm's  master  problem  may  be  effectivelv 
discussed.  Initially,  the  current  estimate  of  the  solution,  f , is 
assumed  to  exist  from  a previous  TR  restricted  master  problem 
solution.  Because  of  the  format  of  the  TR  algorithm,  f may  be 
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expressed  as  a sum  of  N-commodity  flows,  i.e.  f A f^ 


The 


r 1 9 V - 

respective  subproblem  generates  the  set  of  trees,  l<P  ,<p“,  • • j, 
which  minimizes 


Rel  Cost  (f^)  A (T ( f ) , (P-P*)) 


(2.  3.  1) 


for  every  kef  1 , 2, 


Clearly,  Rel  Cost  (pk)  £ 0;  and  if  Rel  Cost  (<P^ ) = 0 for  every 


ke{l,  2,  • • • , nJ  , then  f*  is  not  only  the  restricted  optimal  but  the  global 
optimal. 

The  relative  cost  of  a tree  is  related  to  its  ability  to  reduce  the 
objective,  T(f),  should  it  be  included  in  the  basis.  Assuming  now  that 
the  new  trees  with  negative  relative  costs  have  been  added  to  the  basis. 


the  new  TR  master  problem  requires  a multidimensional  minimization. 

N lk  k 

Specifically,  minimize  T(f)  such  that  f = S L X , <P,.v  where 
- k = 1 i=  1 ik  (l) 


JC  X , = 1.  For  simplicity,  it  is  assumed  that  this  is  done  via  a minor 

1—1  IK 

alteration  on  the  gradient  projection  method  which  Cantor  and  Gerla  use 


for  the  EF  algorithm  [CANT74J.  Calculate  &T(f  )/&Xik  - (?T(f  ),  ( 


for  every  tree  in  the  basis.  Define  Ax  , = -0T(f ' )/QX + (1  /t.  ) • 

« , IK  IKK 

xk 


£ 5T(f  ) /0\  . . Then  it  is  easy  to  show  that  the  direction,  Af, 

m=1  N tk 

defined  by  Af  = £ .£,  A\..  is  a feasible  downhill  direction  from  f 
k = l i = l ik  (i) 


(this  assumes  X-,  > 0 unless  X.,  is  associated  with  a tree  which  has 
ik  ik 


just  been  added  to  the  basis;  see  Appendix  A,  Thm.  2.  3.  A).  The 
condition  > 0 is  maintained  by  dropping  any  tree  from  the  basis  if 
ifs  associated  X^.  is  driven  to  zero  in  the  process  of  the  multi- 


J 


dimensional  minimization. 
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The  next  step  is  to  solve  the  one-dimensional  minimization 

problem:  minimize  T(f  + eAf)  where  e - minimum  (~X  , AX  , )•  For 
0<  e < e:;:  AXik  <0  ik  ik 

the  optimal  value  of  €,  denoted  €,  f f'  + eAf  and  \ , ♦-  x + eA\ ... 

IK.  1 K 1 K 

A new  feasible  downhill  direction  may  be  found  at  this  point,  and  the 
procedure  is  repeated  until  one  is  content  with  the  current  f . 

The  one-dimensional  minimization  requires  about  the  same 
amount  of  effort  whether  it  is  a part  of  the  TR  master  problem  or  of 
the  EF  master  problem.  Moreover,  the  computational  effort  of 

o-  i_ 

finding  Af  is  dominated  by  that  of  finding  (VT(fr),  ) for  every  tree, 

and  then  forming  the  sum,  Af  = - £ ? a>  , <£>!*..  If  t is  the  number 

k=  1 i=  1 ik  (i) 

of  trees  in  the  basis,  then  the  number  of  multiplications  and  additions 
for  the  TR  master  problem  is  about  2t(N-l).  A similar  analysis  for 
finding  Af  in  the  EF  master  problem  results  in  2e(NA)  multiplications 
and  additions  where  e is  the  number  of  extremal  flows  in  the  basis. 

The  efforts  involved  in  finding  Af  seem  comparable,  although  they 
clearly  would  vary  between  particular  problems. 

For  the  number  of  trees  in  the  basis  not  to  grow  out  of  hand,  the 
algorithm  must  remove  from  the  basis  about  as  many  trees  in  the 
master  problem  as  are  added  from  the  previous  subproblem. 

Typically,  only  one  tree  may  have  its  associated  \ , driven  to  zero  for 

1 K 

each  feasible  downhill  direction,  Af,  which  is  found.  Thus,  on  the 
average  only  a small  number  of  trees  with  the  most  negative  costs 
should  be  added  to  the  basis  after  each  subproblem. 

As  it  turns  out  in  the  actual  operation  of  the  algorithm,  the 
potential  problem  of  unbounded  growth  of  the  basis  can  be  avoided.  In 
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the  same  way  that  one  bounds  the  absolute  error  in  any  estimate  of  the 
solution,  one  can  bound  the  reduction  in  the  objective  that  could  be 
achieved  by  an  exact  multidimensional  minimization  with  the  current 
basis.  Specifically,  if 

gk  A minimum  Rel  Cost  (<P,^ .)  1 

““  1 £ i £ t.  {1} 

K 

N k 

then  -Eg  bounds  the  error  due  to  inexact  minimization.  Now  if 

k = 1 

this  bound  is  nearly  as  large  as  the  absolute  bound,  it  is  clear  that  the 
utility  of  adding  columns  is  small  and  more  effort  should  be  spent  on 
the  master  problem  rather  than  in  generating  columns.  From  the 
experience  gained  with  this  algorithm,  it  is  a good  policy  to  add 
columns  only  when  the  ratio  of  the  maximum  error  due  to  inexact 
multidimensional  minimization  to  the  absolute  error  is  under  some 
prespecified  threshold.  Under  this  mode  of  operation,  the  algorithm 
rapidly  reaches  the  point  where  the  rate  of  adding  new  columns  is 
close  to  the  rate  at  which  columns  are  being  evicted. 

In  summary,  the  overall  memory  requirements  and  the  conven- 
ience of  format  of  the  TR  algorithm  would  cause  one  to  favor  it  over 
the  EF  algorithm.  The  relative  speed  of  the  master  problems  is 
difficult  to  predict,  although  it  has  been  shown  that  the  efforts  involved 
in  performing  related  tasks  within  the  master  problems  are  roughly 
equivalent. 

2 . 4 The  CH  Algorithm 

The  CH  algorithm  is  a further  extension  of  the  TR  algorithm  and 

J 
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the  EF  algorithm  in  which  the  decomposition  is  based  on  chains  rather 
than  trees  or  extremal  flows.  The  comparative  discussion  in  this 
section  plus  the  computational  experience  cited  in  Chapter  IV  suggest 
that  the  CII  algorithm  is  competitive  with  respect  to  computation  time 
with  the  EF  and  TR  algorithms,  and  is  competitive  in  terms  of  memory 
requirements  with  the  TR  algorithm. 

This  is  the  master  problem  of  the  CH  algorithm: 

Q 

/N  N rnk  rnl-N 

minimize  T(f)  = T(  L L .Li..  ! 

Vm  = 1 k = l i=l  imk  (0  J 


kAn 


(2.  4.  1) 


such  that 


mk 

L 

i = 1 


imk 


= 1, 


'ik 


* 0 


where  each  <P,.  C is  a member  of  the  given  basis  consisting  of  chains 
(i/ 

(as  defined  in  Section  1.2),  and  c , is  the  number  of  chains  in  the 

mk 

basis  having  origin  m and  destination  k.  In  2.  4. 1,  f is  said  to  be  a 
convex  combination  of  the  chains  in  the  basis.  Again,  denote  the 
solution  of  the  master  problem  as  f . 

The  subproblem  of  the  CH  algorithm  is  analogous  to  the  previous 
subproblems,  namely:  find  the  set  of  chains  which  minimize  (VT(f  ),  f) 
subject  to  1.  1. 3.  If  it  is  found  that  f is  not  the  global  optimum,  then 
the  newly  generated  chains  are  added  to  the  existing  basis  and  the 
algorithm  is  iterated. 

Although  once  again  the  master  problem  has  been  complicated, 
one  hopes  for  a pay-off  due  to  additional  flexibility  in  the  minimization. 
One  may  notice  that  this  algorithm  assumes  the  flavor  of  others  which 
exist  in  the  literature,  in  particular  that  of  Dafermos  [DAFE71].  The 
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principal  difference  between  them  is  in  the  form  of  the  downhill 
directions.  Using  a methodology  similar  to  that  of  the  TR  and  EF 
algorithms  in  the  master  problem,  the  CH  algorithm  finds  a sequence 
of  downhill  directions,  Af,  and  each  generally  alters  the  flow  in  even- 
link  of  the  network;  however,  for  each  Af  in  the  CH  algorithm,  there 
corresponds  N(N-l)  downhill  directions  in  the  Dafermos  algorithm, 
each  of  which  is  associated  with  an  origin-destination  pair  and  a trans- 
fer of  flow  from  the  currently  worst  chain  in  terms  of  marginal  cost 
to  that  of  the  best.  Dafermos'  discussion  of  her  algorithm  presumed 
knowledge  of  every  chain  in  the  network;  however,  clearly  a chain 
generating  scheme  can  be,  in  fact  must  be,  applied  for  networks  of  anv 
size. 


Qualitatively,  the  Dafermos  type  of  scheme  wastes  a great  deal 
of  effort  when  the  algorithm  is  "far"  from  a solution.  One  is  "fine 
tuning"  portions  of  the  network  when  the  flow  in  that  area  may  later  be 
drastically  different  due  to  overall  considerations;  however,  when  the 
entire  network  is  close  to  the  solution,  this  sort  of  sequential  fine 
tuning  may  be  effective.  The  most  likely  meaning  of  the  word  far'  is 
dictated  by  the  point  of  the  occurrence  of  the  tailing  off  phenomenon 
(which  is  well  known  for  the  FW  algorithm  and  also  seems  to  exist  for 
the  EF,  TR  and  CH  algorithms  which  are  no-so-distant  relatives). 
Tailing  off  refers  to  a decrease  in  experimental  convergence  rate  as 
the  algorithm  proceeds.  Although  the  point  of  this  occurrence  is  a 
function  of  many  factors,  such  as  the  demands  on  the  network,  the 
form  of  the  objective  function,  and  the  network  topology,  onset  when  the 
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objective  reaches  from  1 to  2%  of  optimal  is  typical. 

In  terms  of  comparison,  one  can  show  that  a maximum  of 
(NA+N(N- 1 ))  chains  is  required  to  express  the  solution  to  2.  1.  1 (see 
Appendix  A,  Thm.  2.  4.  A).  Chains  have  a maximum  of  (N-l)  nonzero 
terms,  each  of  which  requires  only  the  storage  of  the  arc  number. 

This  means  that  (N-l  )(NA+N(N-  1 ))  bounds  the  memory  requirements 
for  the  basis  of  the  CH  algorithm,  as  compared  to  2(N-1)(NA+N)  for 
the  TR  algorithm.  The  latter  term  is  smaller  in  the  case  that  each 
node,  on  the  average,  has  less  than  (N-2)  nodes  emanating  from  it. 
Clearly,  the  bound  on  the  number  of  nonzero  elements  in  the  chains  is 
quite  liberal.  In  fact,  as  the  network  tends  towards  being  fully 
connected,  the  two  bounds  above  approach  one  another;  and,  since 
chains  in  general  require  even  smaller  storage  areas  than  trees,  the 
CH  algorithm  memory  requirement  becomes  quite  small  compared  to 
the  TR  algorithm.  The  possible  advantages,  however,  are  somewhat 
reduced  by  other  effects  as;  (i)  an  N-fold  increase  in  the  number  of 
identified  commodities  and  the  memory  associated  with  each, 

(ii)  storage  of  variable  length  chains  requires  some  overhead  in 
storage  or  computation  time  since  lengths  are  unpredictable, 

(iii)  storage  is  further  complicated  since  one  desires  the  basis  to  be 
fluid,  i.e.  to  allow  insertion  and  removal  of  chains  at  will,  and  (iv)  the 
ratio  of  the  minimum  number  of  chains  in  a solution  to  the  maximum 
number  of  chains  in  a solution  is  comparatively  large,  so  that  even 
when  the  solution  is  known  to  be  relatively  simple,  a large  percentage 
of  the  maximal  memory  is  still  required. 


I 
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The  extension  of  the  TR  master  problem  to  that  of  the  CH 


master  problem  is  reasonably  straightforward,  and  is  therefore  left 


to  the  reader.  As  in  comparing  the  computational  efforts  of  the  EF 


algorithm  to  the  TR  algorithm  where  analogous  portions  of  the 


algorithms  are  normally  calculated  with  the  same  ease,  the  fundamental 


question  in  comparing  the  TR  algorithm  to  the  CH  algorithm  is  that  of 


the  relative  efforts  of  calculating  Af.  For  the  CH  algorithm,  this 


effort  is  dominated  by:  (i)  calculation  of  (VT(f^),  <P^k)  for  every  chain. 


which  requires  one  multiplication  and  as  many  additions  as  there  are 


nonzero  elements  of<P...  , (ii)  calculation  of  the  AX  , = -ST(f  1/5X-  , 

r ( 1 ) lmk  umk 


+ (l/cm^)  2 ^ 5T(f  )/0X^mk»  which  requires  one  division  for  each 


commodity  and  as  many  additions  as  basis  elements,  and  (iii)  formation 

N N cmk  mk 

of  the  sum  Af  = - 2 2 2 AX.  , <p,..  which  for  anv  given  chain 

k = l m=l  i=  1 imk  (i)  “ g 

requires  one  multiplicalfon  and  as  many  additions  as  nonzero  elements 


of  the  chain.  If  c is  the  number  of  chains  and  n is  the  average  number 


of  links  in  a chain,  then  the  above  sums  to  2c  multiplications  and  N' 

2 

divisions  and  2cn+N  additions.  This  does  not  seem  substantiallv 


different  from  the  2t(N-l)  multiplications  and  additions  of  the  TR 


algorithm. 


2.5  Transitional  Remarks 


A few  remarks  are  in  order  before  concluding  this  chapter.  The 


next  chapter  deals  with  shortest  route  algorithms  which  are  the  tools 
for  solving  the  subproblems  in  an  efficient  manner.  The  effort  associ- 
ated with  solving  the  subproblems  is  nearly  equivalent  to  that  of 
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solving  the  master  problems  [CANT7  -U  , arid  decomposition  techniques 
are  suggested  for  reducing  shortest  route  computation  time  for  certain 
classes  of  networks.  Chapter  IV  discusses  programs  which  solve 
network  problems  with  linear  objectives  which  use  the  same  decom- 
position techniques  of  the  EF,  TR,  and  CH  algorithms,  and  the 
computational  experience  of  that  chapter  will  be  seen  to  support  the 
general  conclusions  of  this  chapter. 
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CHAPTER  III 


SHORTEST  ROUTE  ALGORITHMS 

The  problem  of  finding  all  the  shortest  routes  in  a directed  net- 
work has  an  extensive  literature  [DREY69,  P1ER75]  due  to  the  number 
of  network  problems  to  which  shortest  route  algorithms  are  applied. 
This  chapter  presents  two  new  shortest  route  algorithms  which  can 
significantly  reduce  the  required  computation  time  when  the  network 
is  less  than  fully  connected.  These  algorithms  are  intended  for  use 
where  many  shortest  route  computations  will  be  performed  for  the 
same  topology.  The  first  is  based  on  original  decomposition  ideas 
and  is  called  the  node-by-node  decomposition  (NXN)  algorithm. 

The  second  is  based  on  Hu's  decomposition  algorithm  [HU71,  HU69, 
YEN71]  and  is  designated  the  improved  Hu  (IHU)  algorithm. 

The  shortest  route  problem  is  formulated  as  a shortest  distance 

problem  where  D = [d. . ] is  a given  matrix.  The  number  d. . represents 

3 • 1 

the  length  of  the  directed  arc  from  node  i to  node  j,  and  thus  it  is 

assumed  d. . = 0.  In  this  chapter,  one  further  assumes  that  the  length 

assigned  to  an  arc  is  ® if  and  only  if  the  arc  does  not  exist  in  some 

sense.  A path  P from  i to  j is  an  ordered  sequence  i = kQ,  k^,  • • • , 

km-l'  km  ” j*  and  the  len§th  of  the  path,  L(P),  is  defined  as  L(P)  = 
m 

£ . dk  k . If  P is  any  closed  path,  then  it  is  assumed  L(P)  2 0 
r=l  Kr-lKr 

so  that  the  shortest  distance  problem  is  well  defined.  Then  the 

t After  completion  of  this  work,  the  equivalence  of  the  NXN  algorithm 
with  previous  work  done  at  Network  Analysis  Corporation  under 
ARPA  Order  No.  1523  [NAC71]  was  discovered. 
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problem  is  to  find  D = [d^]  where  d^  min  L(P)  for  P ranging  over 
all  paths  from  i to  j.  Knowing  D alone  does  not  specify  the  shortest 
routes,  but.  it  is  a well  documented  fact  that  by  appropriate  book- 
keeping as  one  calculates  D , the  shortest  routes  can  also  be  estab- 
lished [HU71].  Additio^l  discussion  of  the  problem  of  keeping  track 
of  the  shortest  routes  is  included  in  Section  3.6. 

Typically  D is  calculated  as  a series  of  refinements  on  D. 
Floyd's  algorithm  [FLOY62]  is  cited  for  an  N node  network: 


For  every  ictl, 2, • • • , n) , do  step  a: 

(a)  For  every  j,  ke{l,  2 , • • • , N i , do  step  b; 


(b)  d ..  ♦*  min  (d  ..  , d..  + d ) 
jk  ]k’  ji  lk 


, where  'V  means  "is  replaced  by".  The  algorithm  requires  N‘ 


,3 


additions  and  N'  comparisons,  and  it  is  generally  assumed  additions 


and  comparisons  take  about  the  same  amount  of  time  so  that  one  says 


Floyd's  algorithm  requires  2N  operations.  At  the  conclusion  of  the 
algorithm  D has  replaced  D.  Proof  of  the  algorithm  is  found  else- 
where [HU69],  but  the  interested  reader  can  easily  convince  himself 


that  when  i has  been  stepped  from  1 through  iQ  then  the  current  value 


of  d..  is  the  minimal  distance  over  all  paths  from  j to  k under  the 

jk 


condition  that  the  intermediate  nodes  are  elements  of  the  set 

{l,  2,  • • • , iQ}. 

No  algorithm  which  solves  the  shortest  route  algorithm  could  be 
any  simpler  to  encode,  but  there  are  a variety  of  faster  algorithms  in 
terms  of  number  of  operations  (HOFF76,  SPIR73] . The  standard 
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against  which  the  new  decomposition  algorithms  will  be  measured  is 

Yen's  implementation  of  Dijkstra's  algorithm  [D1JK59,  YEN71] 

3 3 

requiring  g-N  operations.  Those  algorithms  cited  above  claiming 
even  fewer  operations  are  not  significantly  faster,  theoretically,  for 
networks  of  the  size  for  which  computational  experience  is  cited  in 
this  paper;  furthermore,  some  of  the  apparent  gains  of  the 
theoretically  faster  algorithms  would  be  offset  by  their  additional 
algorithmic  complexity. 

3.  1 The  NXN  Algorithm 

The  NXN  algorithm  for  solving  the  shortest  route  problem  is 

3 

actually  a special  case  of  the  following  new  2N  operation  algorithm: 


(1 ) For  every  ie£l,  2,  • • • , N-2}  in  order,  do  step  a: 

(a)  For  every  j,  ke£i+l,  i+2,  • • • , N) , do  step  b: 

(b)  d.,  ♦-  min  (d..  , d..  + d..  ) 
jk  jk*  ji  ik 

(2)  For  every  ic{N-2,  N-3,  • • • , l}  in  order,  uo  step  a: 

(a)  For  every  j,  kefi+1,  i+2,  • * * , n},  do  steps  b and  C; 

(b)  d. . ♦-  min  (d.,  + d,  .,  d.  .) 

i]  ik  kj  ij 

(c)  d..  «-  min  (d..  + d.  d..) 

Ji  jk  ki’  ji 


An  intuitive  proof  of  this  algorithm  will  be  helpful  in  understanding 
the  NXN  algorithm.  By  inductive  reasoning  similar  to  that  for 
Floyd's  algorithm,  when  step  1 has  been  completed  for  i = i , then 
d^  (for  j,  k > iQ)  represents  the  conditional  shortest  j to  k distance 
subject  to  all  intermediate  nodes  being  elements  of  the  set 
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(l,  2,  • • • , i ) . Consequently,  when  step  1 has  been  completed,  then 

the  d.,  (for  j,  k ^ N-2)  represent  unconditional  shortest  distances  d . 
jk  r ij 

Note  that  an  arbitrary  i to  j path  (for  j > i)  must  be  of  the  form 

i.  ' * * . r,  • • * , j , where  r is  the  first  element  in  the  path  such  that 

r > i;  and,  if  this  path  is  the  shortest  path,  then  its  length  is  d.  - d 

r lr  rj 

ij?  **£ 

When  performing  step  2 for  i - N-2,  d^  is  known  and  d.^  must  be  the 

same  as  the  minimal  d conditional  on  all  intermediate  nodes  being 

in  8 

elements  of  the  set  { 1 , 2,  • • • , N-3J;  it  follows  th  t at  the  end  of  step  2 

for  i = N-2,  d.  . = d.  . , and  similarly  d..  = d..  for  everv  ie{N-l,Nj. 
ij  ij  J Ji  Ji  “ J 

Clearly,  inductive  reasoning  shows  that  at  the  end  of  the  algorithm 

D = D . 

The  NXN  algorithm  will  now  be  presented.  However,  in  order 

to  simplify  the  discussion,  it  is  assumed  that  all  of  the  arcs  are 

duplex,  i.  e.  if  d. . < ® then  d..  < ®.  Define  C , called  the  ith 
ij  Ji  i 

connection  set,  as  follows:  jeC.  if  j > i and  there  exists  a path  P from 
i to  j such  that  L(P)  < » and  every  intermediate  node  k satisfies  k < i. 
Notice  that  the  C^  are  functions  of  topology  only. 

In  step  1 of  the  above  algorithm,  d^.  = » if  j<C.  and  d.^  = ® if 
krfC..  Furthermore,  in  step  2 of  the  above  algorithm,  d.,  - ® and 

*■  IK 

^ki  = “ ^kj(C..  The  corresponding  operations  are  clearly  unneces- 
sary; the  algorithm  obtained  by  deleting  them  is  called  the  NXN 
algorithm: 
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(1 ) For  every  ieU,  2,  • • • , N-2  j in  order,  do  step  a: 

(a)  For  every  j,  keC.,  do  step  b: 

(b)  d..  *-  min  (d..  , d..  + d..  ) 

Jk  jk*  ]i  lk 

(2)  For  every  ie{N-2,  N-3,  * * • , 1 } in  order,  do  step  a: 

(a)  For  every  j<[i+lt  i+2,  •••,  n3  and  kcC.,  do  steps 

b and  c: 

(b)  d.  . «-  min  (d..  + d,  .,  d.  .) 

ij  ik  kj*  ij 

(c)  d..  «-  min  (d..  + d,  .,d..) 

Ji  jk  ki*  ]i 

A decomposition  is  defined  as  an  ordering  of  the  nodes.  Since 
the  connection  sets  are  a function  of  the  decomposition,  the  number  of 
operations  which  the  algorithm  requires  is  also  a function  of  the 
decomposition,  as  will  be  demonstrated  in  the  following  section. 

In  the  case  where  some  of  the  arcs  are  not  duplex,  two  alterna- 
tives are  available.  (For  data  communication  networks,  half  duplex 
links  are  rare  because  of  the  resulting  complications  in  error  control.  ) 
The  first  is  to  change  the  definition  of  C.  as  follows:  jeCj  if  j > i and 
there  exists  a path  P from  i to  j or  from  j to  i such  that  L(P)  < * and 
every  intermediate  node  k satisfies  k < i.  This  approach  causes 
unnecessary  operations  for  the  algorithm.  The  alternative  is  to 
define  two  connection  sets  for  each  node- -one  for  the  incoming 
connections  and  one  for  the  outgoing  connections.  In  the  latter  case, 
one  must  alter  the  NXN  algorithm  to  incorporate  the  efficiencies  of 
the  additional  connection  sets.  The  increased  algorithmic  complexity 
of  the  second  approach  and  the  resultant  additional  computer  steps 


f 
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must  be  weighed  against  the  number  of  unnecessary  operations  of  the 
first  approach  for  the  problem  at  bund. 

3 . 2 Decomposing  the  Network  for  the  NXN  Algorithm 

This  section  is  introduced  via  an  example.  Consider  Figs. 

i 

3.2.  1 and  3.2.2  in  which  the  same  network  has  been  decomposed  two 
ways.  For  the  first,  (X  = Ci+1,  N 1,  n}  when  ie(l,  2,  • • • , N-3}  and 
C^-_9  = fN-1,  Nj  ; the  number  of  operations  for  the  NXN  algorithm  is 
calculated  in  a straightforward  fashion  as; 

N-  3 

Step  1,  ' Z (2)  (3)  (3)^  + (2)  (2)  (2) 

v i-  1 / 


/N-3  v 

Step  2,  E (2)  (2)  (3)  (N-i)j  + (2)  (2)  (2)  (2) 

X 1 - 1 / 

9 

which  totals  6N  +12N-66.  By  contrast,  for  the  decomposition  of 

Fig.  3.2.2  C.  - {i+1,  i+2,  • • • , N}  which  is  exactly  the  same  as  if  the 

network  was  fully  connected,  and  it  follows  immediately  that  the  NXN 

3 

algorithm  requires  2N"  operations.  This  example  makes  it  clear  that 
the  choice  of  decomposition  can  have  a profound  effect  on  the  efficiency 
of  the  algorithm. 

For  an  arbitrary  network,  finding  the  optimal  decomposition  in 
the  sense  of  minimizing  the  required  number  of  operations  for  the 
NXN  algorithm  is  not  a trivial  problem  and  probably  can  only  be 
solved  by  exhaustive  comparison.  The  method  of  choosing  the  decom- 
position for  the  examples  which  are  presented  later  in  Section  3.5 
deviated  only  slightly  from  the  following  heuristic  procedure: 
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(1)  Label  a node  "1"  such  that  the  cardinality 
of  is  minimized. 

(2)  For  every  if {2,  3,  • • • , n]  in  order,  do  step  a: 
(a)  Given  the  nodes  which  have  been  labeled 

"1",  ”2",  • • ' , "i-1",  label  an  unlaboled 
node  "i"  such  that  the  cardinality  of  C.  is 
minimized. 


The  effort  in  finding  the  decomposition  via  the  above  procedure  is  on 
the  same  order  as  doing  a shortest  route  computation  via  Floyd's 
algorithm,  and  as  a consequence  computer  time  savings  are  realized 
only  when  the  NXN  algorithm  is  iterated  several  times  for  the  same 

topology. 

There  are  a large  number  of  networks  such  that  the  computation 
time  does  not  vary  widely  with  the  decomposition.  Such  networks 
could  be  termed  "locally  connected"  and  have  the  property  that  the 
nodes  to  which  there  are  direct  arcs  from  any  given  node  are  very 
likely  to  have  direct  arcs  to  one  another.  In  this  case,  the  nodes 
could  be  numbered  very  rapidly  by  eye  with  little  degradation  in 
efficiency  (nodes  must  at  some  point  m time  be  assigned  a number 
anyhow  in  order  to  communicate  the  topology  to  the  computer),  and  in 
the  first  shortest  route  computation  the  connection  sets  could  be 
established  with  very  little  effort.  In  fact,  the  only  modification  to 
the  NXN  algorithm  is  an  additional  step  which  is  included  just  before 
step  la: 
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(aa)  Initially  C.  0;  for  jeti+1,  i+2,  • • • , Nj  , do 
step  (bb): 

(bb)  Include  j in  C.  if  d.  . < ® . 

i ij 

1 i 

The  additional  operations  required  by  this  step  number  N“,  which 
is  quite  modest  for  the  potential  gains. 


3.  3 The  IHU  Algorithm 

The  presentation  of  the  IHU  algorithm  requires  some  additional 

definitions.  For  this  algorithm,  a network  decomposition  is  defined 

as  a division  of  the  network's  nodes  into  ordered  subsets  S, , S„ . • • • . S, 

12  k 

such  that  for  every  i(Sm  and  jeS^,  d„  = =□  if  |m-x|  > 1.  Every  node 

of  the  network  belongs  to  exactly  one  subset.  The  submatrix  D„  <, 

^rm 

contains  all  the  distances  of  arcs  from  elements  of  S.  to  elements  of 
Sm.  and  its  dimensionality  is  |S.  | x |Sm  | (where  Js  | means  the 
cardinality  of  set  S).  Evidently,  „ has  no  finite  entries  in  the 
case  |i-m  | > 1 (see  Fig.  3.3.1). 

Various  matrix  operations  will  be  performed  on  the  submatrices 
to  generate  the  desired  shortest  distance  matrix.  Let  ILL  Q 

bi^i  bjOj 

mean  D„  is  replaced  bv  the  shortest  distance  matrix  computed 
ai°i 

from  the  submatrix  Dc  c . Define  A-B  - (min  (a.  + b .)]  and 

bibi  m ms  mj 

min  (A,  B)  - fmin  (a.  b. .)].  Also  let  S,  U S_  U • • ’ U S = , 

ij  ij  12  mm 

and  if  m = k (where  k is  the  number  of  ordered  sets)  then  Q = 

m 

0^  ^ O.  Define  the  conditional  shortest  distance  submatrix, 

Dc  c.  (0  ),  as  the  shortest  distance  submatrix  under  the  restriction 

OjOj  m 

that  all  the  intermediate  nodes  on  the  respective  conditional  shortest 


Fig.  3.  3.  1 The  form  of  the  D matrix  for  the  IHU 

algorithm  in  the  case  k - 5.  If  the  decom- 
position is  to  be  acceptable,  the  shaded  sub 
matrices  have  no  finite  entries  prior  to  the 
algorithmic  operations  on  the  D matrix. 
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routes  are  members  of  In  the  case  m = k,  „ (Q  ) = 

m S.S.  m 

Ds.s.(n»  4 DS.S.- 
1 J 1 J 

Under  the  assumption  that  an  allowable  decomposition  has  been 
given,  the  following  algorithm  generates  all  the  shortest  distances  in 
the  network  (the  parenthetical  equality  to  the  right  of  each  step  is  tlie 
claim  of  what  each  step  accomplishes); 

(1)  Dc 


S1S1 


4°s  s 

bibi 


(=  d:  (o,)) 

ii  1 


(2)  For  every  ie{l,  2,  •••,  k - 1 } in  order,  do  steps  a,  b, 
c,  and  d; 

(a)  Dc 


'S.  ,S.  DS.  ,S.  ‘ DS.S. 

t+1  i l+l  l ii 


DS  * S.«V> 

l+l  i 


(b)  Ds.s.  Ds.s.  * D 

11+1  11 


'S.S.  (=  DS.S.  /ni)} 

11+1  11+1 


(c)  D 


’s.  ,S.  **  min  (Ds  g , Dg  g • Dg  g ! 

l+l  l+l  l+l  i+I  l+l  i li+l 

(=  D*  q (O.)) 
bi+lSi+l  1 


(d)  Ds.  .S.  *■  *Ds.  ,s 

1+1  1+1  1+1  1+1 


<=D5.+1s.  <°i+l» 

1+1  1+1 


(3)  For  every  ic{k,  k-1,  • • • , 3,  2}  in  order,  do  steps  a, 

t 

b,  and  c; 


(a) 

DS.S. 

*■  Ds.s. 

Ds.s.  <= 

V- 

Ds 

l l-l 

1 1 

1 1-1 

(b) 

Ds.  , s. 

*■  Ds.s. 

* Ds.s.  (= 

* 

Ds 

1-1  1 

1 1-1 

1 1 

(c) 

Ds  s 

bi  i-l 

*-  min  (Dg 

. ,s.*Ds.  ,s. 

1-1  1 1-1  1 

Ds: 

) 


1 1-1 


i-1  i 


<■  Ds.  ,s.  > 

1-1  1-1 


(4)  For  every  re{2,  3,  • • • , k - 1 } in  or  der,  do  step  a: 

(a)  For  every  i,  jc{l , 2,  • • • , k}  if  |i- j | = r,  do  step  b; 


(b)  DS.S.  DS.S  ‘ DS  S. 
1 J i P P .1 


where  p is  an  element  of  the  set  O = {s+  1 , s+2,  • • • , t-2 , t-  1 } for 
s = min(i,  j)  and  t = max  (i,  j)  such  that  |S  | s js^  | for  every  m€Q. 

A rigorous  proof  of  the  algorithm  would  be  very  lengthy  and 
repetitious,  and  the  interested  reader  is  referred  to  Ilu's  work  [HU6P] 
for  exposition  of  a similar  proof.  Steps  1 and  2 are  bootstrapping 
successive  diagonal  and  first  off-diagonal  subn  atrices,  so  that  at  the 

?}e 

end  of  step  2,  Dc  ^ = D . Step  3 is  essentiallv  a backwards  form 

^k^k  ^kbk 

of  step  2 and  replaces  the  diagonal  and  first  off-diagonal  submatrices 
with  the  respective  unconditional  shortest  distance  submatrices. 

Step  4 is  one  method  for  finding  the  unconditional  shortest  distance 
submatrices  corresponding  to  decomposition  sets  which  are  separated 
by  at  least  one  intermediate  set.  The  ordering  in  step  4 allows  p to 
be  any  element  of  the  set  Q,  and  the  particular  choice  of  p minimizes 
the  number  of  operations. 

If  one  assumes  that  the  shortest  distance  calculations  for  sub- 

3 

matrices  are  done  via  Floyd's  method  (requiring  2p  operations  for  a 
p x p submatrix)  and  that  the  pseudo-multiplications  are  done  in  a 
straightforward  manner  (requiring  2pqr  operations  to  calculate  A-B 
where  A is  dimension  p x q and  B is  q x r),  then  the  number  of  opera- 
tions required  by  the  IHU  algorithm  is; 
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Step  1, 

2 |S1  |3 

Step  2, 

(2's1+i!  i 

Is.  1|2  + 

1 1+1  1 

Step  3, 

2 Z (2  |S.  |2  |S.  , | + |S.  . 

i = 2 i i ~ l i-l 

!2  |s.  1) 

i ■ 

Step  4,  2 £ |S.  I |S.  I |S  I 

i,  j i J P 

|i“  j I>1 


The  total  number  of  operations  is  then: 


2 (kL  |S.  u S.  |3 
\ 1=1  i l+l 


k- 1 

& |s|‘ 


+ E 
i,  j 

H >i 


Is. I |Sp 


One  may  compare  the  JHl1  algorithm  to  other  versions  of  Hu's 
algorithm.  For  any  given  decomposition,  the  IHU  algorithm  requires 
fewer  operations  than  the  fastest  version  of  Hu's  algorithm  known  to 
the  author,  which  is  that  due  to  Yen  [YEN71] . For  purposes  of 
comparison,  an  example  which  commonly  appears  in  the  literature 
[HI1 71,  HUB 9,  YEN7 1 ] is  presented.  Let  |S.  | = $ for  i even  ana  |S.  | = t 
for  i odd.  Assume  ^ X,  and  let  k,  the  number  of  sets,  be  odd. 

Define  m - (k+l)/2.  In  this  case,  the  new  algorithm  requires 

2(mt  + (m  + 5m  - Blt^g  + (2m“  + 2m  - 6)t5~  + (m2  - 4m  + 5)5  ?)  operations 

Yen's  modification  requires  2(mt3  + (m2  + 6m  - 7)t25  + (2m2  + 10m  - 20)t5 

2 3 

+ (m  + 6m  - 14)5  ).  The  new  algorithm  is  faster  for  the  entire  range 
of  interest,  i.e.  t 2 5 2 1 and  m 2 2.  As  a particular  case,  let  5 = t 
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i 


and  m =3;  the  IHU  algorithm  requires  82t ^ operations.  Yen's 

modification  requires  128t  operations,  and  Floyd's  algorithm  reqv 
3 

250t  operations. 

3.4  Decomposing  the  Network  for  the  IHU  Algor ithm 

Perhaps  even  more  important  than  the  numerical  gains  of  the  n.  > 
algorithm  are  the  insights  it  provides  into  optimal  decomposition  of  a 
network.  Assume  that  Floyd's  method  is  used  for  shortest  route 

M 

computations  on  submatrices,  and  that  pseudo-multiplications  are 
done  by  the  straightforward  technique.  It  follows  that  for  a given 
decomposition,  if  a further  decomposition  exists  by  partitioning  of 
existing  sets,  then  the  computation  time  of  the  further  decomposition 
is  less  than  that  of  the  given  decomposition.  This  "more  the  better" 
fact  suggests  a heuristically  good  decomposition  technique  which  can 
be  performed  by  the  computer  or  quickly  guessed  at  by  eye.  If  the 
decomposition  is  to  be  done  automatically  by  the  computer,  however, 
it  should  probably  be  limited  to  those  cases  where  many  shortest 
route  computations  for  the  same  topology  will  be  performed,  as  in 
column-generating  linear  programs.  An  algorithm  for  finding  a good 
network  decomposition  for  the  IHU  algorithm  is; 

(a)  find  two  nodes,  j and  k,  such  that  the  minimal 
number  of  arcs,  d,  connecting  them  is  maximal 
over  all  pairs  of  nodes;  i.e.  find  the  diameter 
of  the  network  and  an  associated  pair  of  nodes; 
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(b;  construct  d+1  sets  by  letting  { jj  and 

S.+1  = {m  Imefo-flJ  and  drm  < ® for  some  rcSji  . 

This  procedure  was  used  to  generate  the  IHU  decomposition  sets  for 
the  examples  of  the  next  section,  and  the  reader  may  want  to  look  at 
the  figures  associated  with  that  section  at  this  point. 

3.  5 Some  Examples  Using  the  IHU  and  NXN  Algorithms 

In  this  section,  several  examples  are  presented  which  provide 
insight  into  the  classes  of  networks  for  which  the  NXN  and  IHU 
algorithms  can  substantially  reduce  shortest  route  computation  time. 
One  can  probably  deduce  from  these  examples  that  the  NXN  algorithm 
will  perform  a shortest  route  computation  more  rapidly  on  a "typical" 
network.  However,  the  fourth  example  suggests  that  there  is  a small 
and  special  class  of  networks  for  which  the  IHU  algorithm  is  more 
efficient.  Typically,  such  networks  may  be  decomposed  in  such  a 
manner  as  to  be  a variation  on  the  following  theme:  |S  | for  i odd  is 

large  compared  to  |s.  | for  i even,  and  if  jeS.  and  k«S.  then  j and  k are 
very  likely  to  have  direct  arcs  to  one  another. 

The  first  example  is  an  old  version  of  the  ARPA  net  which  is 
shown  in  Fig.  3.5.  1.  In  that  figure,  the  NXN  decomposition  is 
defined  by  the  numbering  of  the  nodes,  and  the  IHU  decomposition  is 
defined  by  the  partitioning  of  the  nodes  with  broken  lines.  This  net- 
work lends  itself  to  NXN  decomposition  due  to  the  large  number  of 
nodes  which  have  arcs  directly  to  only  two  other  nodes--a  fact  which 
keeps  the  cardinality  of  connection  sets  very  low. 


The  second  example  is  the  17-node  symmetric  network  shown 
in  Fig.  3.5.2,  This  network  is  not  "locally  connected"  to  a very 
large  degree,  but  still  the  NXN  algorithm  is  (perhaps  surprisingly) 
efficient. 

The  third  example  is  the  64-node  network  displayed  in  Fig.  3.5.  3. 
The  density  of  arcs  is  perhaps  greater  here  than  in  the  other  examples, 
but  a large  degree  of  local  connectivity  promotes  the  efficiency  of  the 
NXN  algorithm. 

The  fourth  example  is  the  36-node  network  of  Fig.  3.5.4.  This 
is  typical  of  the  type  of  network  for  which  one  would  expect  the  IHU 
algorithm  to  be  more  efficient  than  the  NXN  algorithm. 

Efficiency  is  measured  with  Yen's  implementation  of  Dijkstra's 
algorithm  as  the  standard.  Theoretical  efficiency  refers  to  the 
relative  savings  in  the  number  of  operations  required  to  perform  a 
shortest  route  computation.  The  computation  times  for  the  IBM 
370-168  to  execute  the  Fortran  program  of  the  various  algorithms 
were  noted,  and  relative  savings  are  referred  to  as  the  measured 
efficiency.  The  comparisons  of  the  various  algorithms  in  performing 
shortest  route  calculations  on  the  three  sample  networks  are  summar- 
ized in  Table  3.  5.  1.  The  Fortran  programs  were  compiled  by  the 
IBM  G1  compiler;  and  each  algorithm  not  only  computed  the  shortest 
distance  matrix,  but  also  computed  a routing  matrix  which  specified 
the  first  arc  from  each  node  on  the  shortest  route  to  any  other  node. 


Fig.  3.5.2  A 47  node  symmetric  network  (nodes  are  connected  to  first  and 
seventh  nearest  neighbors  by  arcs).  XXN  decomposition  is 
indicated  by  node  labeling..  IHU  decomposition  sets:  Si  = { 1 }, 
S2  = 121.47,28,  39},  S3  = 134,  40,  30,  26,  2,  4,  5,  8],  S4 
(14,  17,  18.29,22,23,  42,  41, 32,  25,  36,  43],  S5  = (46,  33,  15,  37, 
31,  38  44,  3,  6,  7,  9,  29,  1 1.  12},  and  Sfi  = {l  3.  27,  16,  19,  IS,  35, 
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Fig.  3.  5.  4 A 36  node  network.  NXN  decomposition  is  indicated 
by  node  labeling.  IHU  decomposition  sets; 

Sj  = (1,2,3,  4,  5,  6.  7,8,9}  , S2  --  £l0,  11  j , S3  = {l2,  13, 
14,  15,  16,  17,  18},  S4  - {19,20},  S5  = {21,  22,2  3,  24, 

25,  26,  27} , Sg  = {28,  29}  , S?  = {30,  31,  32,  33,  34,  35,  36; 


ARPA 
network 
of  Fig. 
3.  5.  1 

47 -node 
network 
of  Fig. 
3.5.2 

64-node 
network 
of  Fig. 
3.5.3 

36 -node 
network 
of  Fig. 
3.  5.  4 

Dijkstra's 

shortest 

route 

algorithm 

Number  of 
operations 

27040 

154630 

393216 

69984 

Theoretical 

efficiency 

1. 00 

1. 00 

1. 00 

1. 00 

Computation 
time  (seconds) 

0.  090 

0.  450 

1.115 

0.210 

Measured 

efficiency 

1. 00 

1.00 

1.  on 

1. 00 

IHU 

algorithm 

Number  of 
operations 

6948 

83880 

124722 

11360 

Theoretical 

efficiency 

3.  89 

1.84 

3.  15 

6.  16 

Computation 
time  (seconds) 

0.  020 

0.  195 

0.290 

• 0.  040 

Measured 

efficiency 

4.  50 

2.31 

3.84 

5.25 

NXN 

algorithm 

Number  of 
operations 

2828 

28608 

42416 

13316 

Theoretical 

efficiency 

9.  56 

5.41 

i 

9.27 

5.26 

Computation 
time  (seconds) 

0.  015 

0.  1 15 

0.  165 

0.  050 

Measured 

efficiency 

6.  00 

3.  91 

6.  76 

4.20 

Table  3.5.1  Comparative  performance  of  three  different  shortest 
route  algorithms  on  the  three  sample  networks. 
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3.  6 Additional  Notes  on  the  Shortest  Route  Algorithms 

For  completeness,  the  method  of  bookkeeping  for  finding  the 

shortest  routes  which  was  implemented  in  the  algorithms  encoded  as 

a part  of  this  thesis  is  presented.  An  N x N "next  arc"  matrix, 

X = fx  . .1,  ni  established,  where  x . is  the  first  arc  on  the  current 
hi  ‘j 

estimate  of  the  shortest  route  from  i to  j.  Initially,  one  sets 
x_  = m where  OR(m)  = i and  DS(m)  = j for  every  existing  arc  m. 
Whenever  an  operation  d.j  ♦-  min  (d.j,  + d^  ) is  performed  such  that 

d.^  + d^  is  the  distinct  minimum,  then  one  makes  the  replacement 
x..  «-  x..  . At  the  conclusion  of  a shortest  distance  algorithm,  the 
shortest  routes  may  be  generated  in  an  obvious  fashion. 

As  a practical  matter,  it  is  noted  that  there  are  cases  where  a 
shortest  distance  algorithm  succeeds  in  finding  the  shortest  distances, 
but  the  bookkeeping  for  shortest  routes  fails.  Failure  is  even  more 
likely  in  a digital  computer  implementation  than  in  an  exact  calcula- 
tion (although  both  are  possible).  Sufficient  conditions  for  success 
of  the  bookkeeping  are  as  follows:  (i)  the  network  is  connected,  i.  e. 

there  exists  a finite  distance  path  between  every  ordered  pair  of 
nodes,  and  (ii)  if  d is  the  minimum  distance  assigned  to  any 

existing  arc  and  d is  the  maximum  distance  assigned  to  any 

6 max  6 

existing  arc,  then  d > 0 and  d + d is  distinct  from  d in 

min  min  max  max 

the  computer  storage.  Given  these  conditions  as  a hint,  the 
interested  reader  should  be  able  to  construct  an  example  where  the 
bookkeeping  fails  although  a shortest  distance  algorithm  succeeds. 

The  use  of  the  NXN  and  IHU  algorithms  should  be  limited  to 
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networks  of  moderate  size,  i.e.  under  several  hundred  nodes.  The 
basic  reason  for  this  limitation  is  that  there  are  several  algorithms 
that  do  not  require  storage  of  the  entire  D matrix  and  are 
consequently  more  efficient  for  very  large  sparse  networks.  Two  of 
these,  Dijkstra's  algorithm  (not  Yen's  version)  and  Bellman’s 
algorithm,  were  compared  by  Golden  [GOL.D76].  From  Golden’s 
work,  one  can  see  that  for  large  sparse  networks  the  algorithms  he 
compared  are  quite  efficient;  their  performance,  however,  is  not 
nearly  so  impressive  for  smaller  networks.  The  failure  of  these 
algorithms  which  do  not  explicitly  manipulate  the  D matrix  to  perform 
competitively  on  the  smaller  networks  is  owed  to  a large  number  of 
overhead  operations  associated  with  such  algorithms. 


CHAPTER  IV 


A LINEAR  MULTICOMMODITY  ROUTING  PROBLEM 

This  chapter  studies  some  methods  for  solving  a linear  multi- 
commodity  routing  problem  which  yield  a solution  that  has  certain 
desirable  properties  in  the  context  of  computer  networks.  In  essence, 
the  solution  to  the  linear  routing  problem  places  the  minimum  possible 
demand  on  the  channels  which  form  the  "bottleneck"  of  the  network. 

The  linear  problem  is  to  find 

min  ( max  f.  / C . ) 
fcF  1 si  SNA  1 1 

That  is,  one  is  minimizing  the  maximum  utilization  factor  of  any 
channel  in  the  network  over  all  possible  routing  strategies. 

The  channels  associated  with  the  binding  capacity  constraints  in 
the  solution  to  the  linear  program  form  the  bottleneck  of  the  network. 
Binding  constraints,  in  this  chapter,  are  those  which  are  necessarily 
tight  or  active  at  the  solution  of  the  problem.  They  are  "binding"  in 
the  sense  that  a decrease  of  the  capacity  of  any  associated  channel 
implies  an  increase  in  the  value  of  the  objective.  There  are  no  flows 
through  the  bottleneck  in  the  solution  which  even  in  part  could  avoid 
the  bottleneck,  and  it  is  in  this  sense  that  the  demand  on  the  bottleneck 
is  minimized. 

From  a slightly  different  point  of  view,  this  solution  solves  a 
throughput  problem  for  the  network.  The  inverse  of  the  value  of  this 
linear  objective  function  at  the  solution  is  the  maximum  scale  factor 


k. 
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on  the  entire  demand  matrix  such  that  a feasible  flow  (one  which  does 
not  exceed  any  channel  capacity)  exists. 


4.  1 Decomposition  Techniques  for  Solution  to 
the  Linear  Routing  Problem 

The  outlook  of  this  chapter  will  be  from  the  viewpoint  of  the 
Dantzig-  Wolfe  decomposition  technique  [DANT63].  The  same 
approach  could  have  been  taken  as  in  Chapter  II,  i.e.  viewing  the  tech- 
niques as  column  generating  algorithms.  The  problem  which  is  to  be 
solved, 

f. 

(4.1.1) 


min  ( max  -4- ) 
f€F  1 ^NA  i 


may  be  restated  as: 


minimize  such  that 

-^C  + If1  + If2  + ir 


If1'1  * 0 


Ef 


= bJ 


Ef 


= b‘ 


(4. 1.2) 


EC 


= b‘ 


Ef3^  = bN 


* 0 for  kc{l,  2.  • • *,  NJ 


where  C = [Cj,  Cg,  • • • , is  a real  scalar,  I is  the  NA  x NA 

identity  matrix,  and  the  remainder  of  the  notation  is  consistent  with 


1. 1.2. 

According  to  the  Dantzig- Wolfe  decomposition  principle,  one 
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may  formulate  the  master  problem  of  4.  1.2  as: 


minimize  (3  such  that 

-^1C  + + JiXj2^j) 


+ f.l  5 ° 


tl 

£ x 

j=i  Jl 


-2 

£ X >9 
j = l J2 


= 1 


N 


£ X_  = l 


3=1 


']N 


Xjk  s 0 for  every  j € C 1 , 2, 


. t^i  and  kcfl,  2, • • • , n) 


(4.  1.3) 


where  the  are  solutions  to  the  various  subproblems  which  will  be 
defined  momentarily,  and  t,  is  the  number  of  trees  associated  with 
the  destination  node  k.  The  above  format  should  be  compared  with 


that  of  2.2.1.  Upon  solution  to  the  restricted  master,  4.1.3,  there 


exists  a set  of  dual  variables,  ' -J  -j  = • • * • ffNA  ’ ffl  * CT<>  ’ * • • »aNJ' 


associated  with  the  constraints  of  the  master;  77.  is  the  dual  variable 


associated  with  the  ith  arc,  and  cr  is  the  dual  variable  associated  with 
the  constraint  £ \ = l. 

j = l Jk 


There  are  N subproblems  which  are  solved  after  each  master, 
each  having  the  form. 


J 


-c; 


minimize 


-7 V ^ 


such  that 


Ef^  = bk,  fk  ^ 0 


(4.1.4) 


and  the  solution  is  denoted^  + j y Each  + ^ is  a column  which 

is  added  to  the  current  master  if  its  cost  relative  to  the  current 

solution  is  negative,  i.e.  if  (-TT'cr^  + ^ - a^)  < 0.  In  the  case  that  the 

k 

column  is  added,  t,  *-  t + 1 in  the  master  problem.  It  is  noted  that 
4.1.4  is  the  same  problem  as  1.2.2,  so  that  the  method  of  solution  is 
via  shortest  route  algorithms.  Consequently,  the  solutions  to  4.  1. 4 
are  trees.  It  is  known  from  linear  programming  that  the  dual 
variables  associated  with  the  arcs  are  nonpositive,  assuring  the 
success  of  shortest  route  algorithms  for  solving  the  subproblem 
[DANT63].  If  all  of  the  N trees  generated  by  the  subproblems  have 
nonnegative  relative  cost,  then  the  solution  to  4.  1.2  was  established 
by  the  previous  master  problem  and  the  algorithm  terminates. 
Otherwise  the  master  problem  is  resolved,  and  the  above  procedure  is 
repeated. 

The  number  of  constraints  in  4.  1. 3 is  NA+N,  which  is  also  the 
maximal  number  of  trees  in  a solution.  the  value  of  fi  at  the 

solution  of  4.1.3,  is  the  minimal  scale  factor  for  the  capacities  such 
that  there  exists  a feasible  routing  strategy,  i.e.  if  6 < $ then  there 
does  not  exist  an  feF  such  that  f.  s gc.  for  every  ict  1 , 2,  • • • , NAj. 

The  above  algorithm  is  denoted  the  LTR  algorithm  (linear 


i 


routing  algorithm  based  on  a tree  decomposition).  Alternatively, 
4.  1.2  could  have  been  expressed  as: 


r 


-tib- 


minimize  (i.  such  that 


-/3  C + If  s 0 
Mf  ^ b 
f ^ 0 


(4.  1.5) 


where  the  notation  is  consistent  with  1.1.1.  This  problem  may  be 
solved  via  the  Dantzig -Wolfe  decomposit  on  principle,  and  the 
procedure  is  similar  to  that  of  the  EF  algorithm  which  was  presented 
in  Section  2.  1.  The  procedure  is  called  the  LEF  algorithm,  and  the 
master  problem  is; 

minimize  /3  such  that 


-hc  s 0 


I 1,  ■ 1 

i = 1 1 


\ 20 


(4.  1.6) 


with  a subproblem  analogous  to  that  of  4.  1.4  and  a similar  criterion 
for  terminating  the  algorithm.  The  solution  to  4.  1.2  has  no  more 
than  NA+1  extremal  flows  in  the  solution. 

9 

Finally,  4.  1.2  may  be  formulated  as  an  N“-commodity  flow 
problem  using  the  notation  of  1.  1.3.  Similarly,  a decomposition 
scheme  for  solving  the  problem  based  on  chains  can  be  derived,  and 
the  resultant  algorithm  is  called  the  LCII  algorithm.  Accordingly, 
the  solution  to  4.1.2  would  have  no  more  than  NA+N(N-1)  chains  in  the 
solution.  Note  that  upon  finding  the  solution  to  4.  1.  1,  the  associated 
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routing  strategy  is  derived  in  the  same  manner  as  presented  in 
Section  2.  1 with  varying  amounts  of  difficulty  according  to  the  decom- 
position scheme  and  its  format  for  the  solution. 

4.  2 Higher  Order  Criteria  for  the  Linear  Program 
Routing  Strategy 

This  section  discusses  higher  order  criteria  for  the  solution  to 
4.1.1.  There  are  many  members  of  F which  solve  4.  1.  1,  and  the 
purpose  of  the  criteria  is  to  define  the  "best"  among  the  many  solutions. 
For  the  sake  of  compact  notation,  the  criteria  are  exemplified  in  terms 
of  the  LEF  algorithm,  with  the  extensions  to  the  LTR  algorithm  and 
the  LCH  algorithm  being  sufficiently  straightforward  as  to  be  omitted. 
Assume  that  4.  1. 1 has  been  solved  by  the  LEF  algorithm  and  its 
solution  is  p . 

Having  solved  4.  1. 6,  some  of  the  tight  capacity  constraints  are 

associated  with  a dual  variable  which  is  negative.  It  is  an  immediate 

consequence  of  duality  [DANT63,  pp.  134-40]  that  such  a constraint 

has  no  slack  in  any  solution  of  4.  1. 6;  hence,  any  reduction  of  the 

capacity  of  the  associated  channel  would  result  in  an  increase  in  the 

objective  function.  Of  course,  every  channel  in  the  bottleneck 

(i.e.  every  channel  corresponding  to  a binding  constraint)  is  not 

necessarily  indicated  by  a negative  dual  variable;  however,  provisions 

for  such  omissions  are  included  in  the  following  discussion. 

Define  (called  the  first  saturation  set)  as  follows;  itS^  if 

ff.  < 0 in  the  solution  of  4.1.6.  Also,  define  g.  Q = 1 if  icS  , and 
1 i#  1 

6-  o = 0 if  ifS  . Then  the  secondary  criterion  for  the  solution  to  the 

1,  Oj  l 
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linear  program  ' s defined  to  be  that  the  solution  to  4.  1 . 1 also  sol 

minimize  j39  such  that 
+ If  s ^C1 
Mf  ^ b 


f ^ 0 

where  C1  = [c,  ..c,  „,  • • • , c A ] ' such  that  c.  . = C.5.  e , and 

1,1  1,  z 1,  NA  1,1  l i,bj 

where  C1  = [c^  v 2,  • • • . NA] ' such  that  c^.  - Cjd-fi^). 

Note  that  C*  + C*  = C.  The  physical  interpretation  of  the  problem  is 
to  minimize  the  maximal  utilization  factor  of  any  channel  other  than 
those  in  given  that  the  utilization  factor  of  the  channels  of  S ^ are 
fixed  at  p^.  If  the  value  of  the  objective  of  the  secondary  problem  is 
equal  to  p , then  the  entire  bottleneck  is  not  indicated  by  S^.  In  this 
case  the  channels  associated  with  negative  dual  variables  in  the 
solution  to  the  secondary  problem  are  also  part  of  the  bottleneck. 
However,  computational  experience  indicates  that  if  the  elements  of 
the  demand  matrix  are  arbitrarily  chosen  and  few  are  zero,  then  it  is 
fairly  likely  that  indicates  the  entire  bottleneck. 

The  jth  order  criterion  is  a straightforward  generalization  of 
the  secondary  criterion.  The  reader  will  be  spared  the  additional 
notational  complexities,  although  the  interpretation  is  as  follows. 

For  every  criterion,  there  are  constraints  in  the  solution  of  the 
associated  problem  which  are  for  the  first  time  assigned  negative  dual 
variables.  The  utilization  factors  of  the  corresponding  channels  are 
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then  fixed  at  their  value  in  the  solution.  The  problem  associated  with 
the  next  criterion  is  to  minimize  the  maximal  utilization  factor  over 
those  channels  which  have  not  had  their  utilization  factors  fixed  bv  a 
previous  criterion. 

4.  3 Examples  Using  the  LEF,  LTR  and  LCH  Algorithms 

Several  examples  which  provide  insight  into  the  behavior  of  the 
linear  programming  algorithms  are  presented  in  this  section.  A 
linear  programming  system  entitled  SEXOP  (Subroutines  for  Experi- 
mental Optimization),  which  is  currently  maintained  at  MIT  by 
Prof.  Roy  Marsten  [MARS74],  was  employed  for  these  examples  in 
the  master  problem.  Some  limitations  on  the  size  of  the  examples 
were  imposed  by  the  capabilities  of  SEXOP.  The  principal  constraints 
were  that  no  network  with  more  than  99  directed  channels  could  be 
handled,  and  the  LCH  algorithm  could  not  be  used  on  networks  with 
more  than  14  nodes. 

The  four  example  networks  are  shown  in  Figs.  4.3.1,  4.3.2. 
4.3.3,  and  4.3.4.  The  requirement  matrix  for  each  of  the  problems 
was  pseudo- randomly  generated  according  to  the  rules  of  Appendix  C. 
The  arcs  are  assumed  to  be  duplex  and  have  unit  capacity.  The  arcs 
marked  by  arrows  in  the  figures  denote  the  channels  associated  with 
the  binding  constraints  in  the  solution  of  the  specified  problem.  The 

A 

algorithms  successively  refine  the  estimate  of^,  i.e.  0 until  a 
solution  is  obtained,  and  graphs  of  0.  versus  computation  time  are 
presented  in  Figs.  4.3.5,  4.3.6,  4.3.7,  and  4.3.8. 


V 


,1  .2.  .3  A-  5 .to  T S 9 IO  seconds 

Fig.  4.3.5  Experimental  convergence  of  of  the  linear 
routing  problem  associated  with  Fig.  4.3.  1 
for  the  LEF,  LTR  and  LCH  algorithms. 
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Fig.  4.3.6  Experimental  convergence  of  0 of  the  linear 
routing  problem  associated  with  Fig.  4.3.2 
for  the  LEF  and  LTH  algorithms  (the  LCH 
algorithm  was  not  run  for  this  problem). 
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Fig.  4.3.7  Experimental  convergence  of  /9  of  the  linear 
routing  problem  associated  with  Fig„  4.3.3 
for  the  LEF,  LTR  and  LCH  algorithms  (note 
the  change  of  horizontal  scale  at  3 seconds). 


routing  problem  associated  with  Fig.  4.3.4  for 
the  LEF,  LTR  and  LCH  algorithms.  The  LEF 
algorithm  still  had  not  reached  the  solution  at 
90  seconds  when  its  value  of  was  6.  00. 
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It  is  noted  that  in  the  case  of  unit  capacity  for  every  arc,  the 

linear  routing  problem  becomes;  min  max  f . Thinking  of 

fcF  1 si  SNA  1 

channels  in  the  queueing  sense,  and  the  problem  as  minimizing  the 
maximal  utilization  factor,  then  values  of  the  objective  under  one  are 
desired.  This  is  not  the  case  for  these  problems,  as  is  displaced  in 
the  graphs,  but  one  may  scale  the  capacities  in  order  to  achieve 
feasibility.  Furthermore,  the  graphs  only  display  the  convergence 
for  the  first  criterion.  Refining  the  solution  for  higher  order 
criteria  was  found  to  consume  large  amounts  of  computation  time 
relative  to  that  of  the  first  order  criterion  (see  Appendix  D). 

The  comparative  speeds  at  which  the  various  algorithms  con- 
verge to  the  solution  can  be  explained  qualitatively  as  follows.  The 
principal  component  in  the  computation  time  of  a linear  program  is 
that  associated  with  the  construction  and  maintenance  of  the  basis 
inverse.  Each  time  a column  is  added  to  the  basis,  the  basis  inverse 
must  be  altered.  Naturally,  one  would  like  to  minimize  the  number 
of  changes  in  the  basis.  Consider  being  very  "far"  from  the  solution. 
The  addition  oi  one  chain  to  the  basis  will  probably  have  a very  modest 
effect  on  the  objective.  The  addition  of  a tree  would  have  more 
power,  and  an  extremal  flow  even  more.  However,  as  the  algorithm 
closes  in  towards  the  solution,  the  flexibility  of  the  columns  becomes 
the  most  important  factor  in  the  convergence,  and  the  three  types  of 
columns  would  be  desirable  in  the  reverse  order  of  'heir  power  during 
the  initial  stages  of  the  algorithm. 

The  conclusion  is  that  one  would  expect  the  LEF  algorithm  to 
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be  the  fastest  starter  and  the  LCH  algorithm  to  be  the  fastest  finisher. 
If  the  problem  is  very  complicated  in  the  sense  that  many  constraints 
are  binding  in  the  solution,  then  the  finish  is  far  more  important  than 
the  start  and  the  LCH  algorithm  will  be  the  clear  winner.  If  the 
problem  is  not  very  comnlicated,  then  a faster  starting  algorithm  mav 
also  be  the  first  one  to  solve  the  problem  (see  Fig.  4.  3.  5). 

The  conclusions  from  the  linear  programming  examples  are 
fairly  clear.  The  LEF  algorithm  (which  is  very  much  related  to  the 
method  Hu  suggests  for  solving  the  problem  [HU69])  will  rarely  be  the 
best  technique  for  solving  the  problem.  The  LTR  algorithm  has  a 
compact  representation  in  the  machine  and  has  good  speed  character- 
istics relative  to  the  LEF  algorithm  and  should  nearly  always  be 
preferred  to  the  LEF  algorithm.  The  only  drawback  of  the  LCH 
algorithm,  which  is  fastest  for  the  more  difficult  problems  where 
speed  is  most  important,  is  the  very  large  number  of  commodities. 
The  columns  of  the  LEF  and  LTR  algorithms  have  a uniform  number 
of  nonzero  elements.  A general  purpose  linear  programming  package 
like  SEXOP  handles  the  storage  of  the  columns  efficiently.  However, 
the  unpredictability  of  the  number  of  nonzero  elements  of  the  columns 
of  the  LCH  algorithm  compounded  by  the  failure  to  exploit  the  fact  that 
all  the  nonzero  elements  in  any  column  are  equal,  makes  storage  in 
the  SEXOP  package  wasteful.  A special  purpose  linear  program 
would  alleviate  some  of  the  storage  burdens. 
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4.  4 A Nonlinear  Programming  Approach  to  the 

Linear  Routing  Problem 

This  section  opens  with  a presentation  of  an  algorithm  which  is 
based  on  nonlinear  optimization  techniques,  but  which  converges  to 
the  solution  on  the  linear  routing  problem.  As  an  aid,  a lower  bound 
on  ^ is  derived  which  may  be  used  to  make  an  intelligent  decision  on 
when  to  terminate  the  algorithm. 

Prior  to  a formal  statement  of  the  algorithm,  the  intuition 
behind  the  algorithm  is  presented.  The  parameter  a is  thought  of  as 
a scaling  factor  on  the  capacity  vector  of  the  network.  An  arbitrary 
fcF  is  selected  as  the  starting  point  for  the  algorithm.  An  initial  cr 
is  chosen  such  that  f is  feasible,  i.  e.  f <crC.  A minimization  is 
performed  with  an  objective  function  which  has  a nole  at  each  f = aC.. 
This  tends  to  drive  down  the  higher  values  of  f.  / C.  at  the  expense  of 
the  others.  Then  a is  rechosen  such  that  the  slack  in  the  problem  is 
reduced.  However,  the  new  value  ofe  must  not  destroy  the  feasibility 
of  the  current  flow  vector  (which  is  thought  of  as  the  current  estimate 
of  the  solution  to  4.  1.  1).  The  algorithm  either  stops  according  to 
some  criterion  or  iterates  by  returning  to  the  minimization  step. 

The  algorithm.  The  following  algorithm  will  be  denoted 
AL*(a,  e,y): 

(1)  Choose  arbitrarily  any  f^cF.  Define  /3  (f)  = max  (f./C.), 

Q 1 1 si  SNA  1 1 

and  set  j(f  ) where  y has  been  chosen  greater  than  1. 

Let  K be  an  iteration  counter,  and  set  K = 1. 


r 


: 

i 


(2)  To  within  a predetermined  tolerance,  solve  the  following 


-81- 


problem  (starting  from  the  flow  vector 


-1, 


min  P(u,.,  f) 
feF  K 


where 


P(cr,  f) 


NA  <*C. 


r 

u 

j = i 


'3  ’ *3 


NA 

Li 

j=i 


TT75- 


and  where  the  utilization  factor  f./C^  is  denoted  by  p..  That 
is,  find  an  such  that 

P(<T  f1^)  s (1  + 0 min  P(or  , f) 

K f€F  K 

where  c has  been  chosen  as  an  arbitrary  positive  constant. 
This  tolerance  may  be  guaranteed  by  the  upper  bound  on  the 
error  given  in  Chapter  II. 

(3)  Set , = (l-a)CT..  + aft  (f^),  assuming  a has  been  chosen 
K+ 1 K 1 

such  that  0 < a < 1.  If  K = K , where  K is  a prespecified 

o o 

number  of  iterations,  cease  the  algorithm.  Otherwise, 
perform  K «-  K + 1,  and  go  to  step  2. 

From  Chapter  II,  one  has  a choice  of  several  algorithms  to  use 
in  step  2.  Each  performs  the  required  step  in  a finite  number  of 
operations.  If  the  EF  algorithm  is  employed,  the  above  algorithm  is 
denoted  EF^a,  c ,y),  and  analogous  notation  is  used  for  the  other 
algorithms. 

Because  of  the  finite  effort  involved  in  step  2,  it  is  now  shown 


that  AL*(a,  e,y)  converges  to  a 6 tolerance  solution  of  the  linear 
routing  problem  in  a finite  number  of  operations;  i.e.  for  any  0, 


there  exists  a Kq  < =o  such  that  ^(f^0)  s 5 + ,3^.  The  proof  is  not 
difficult.  Consider  step  2.  Let  g be  any  solution  of  the  linear 
routing  problem.  Then 


a NA 

min  P(a  f)  s P(a  g)  x ■ 

feF  K a..  -B, 

tx  ■ 1 


since  each  element  of  g satisfies 

P(tfK,fK)  * 


gj/Cj  = 

NA(  1 + e)aK 


~K 


This  implies 


Clearly.  V Hx  ~ ^1  (fK))  £ P(aK>  fK)  since  the  former  is  one  term  of 
the  summation  which  defines  the  latter,  and  it  follows  that 

a NA(1  + €)cr 

ix s K 

c-K-^(fK)  ' Tk-S, 

Note  that  cr^.  > ^(f^)  s ^ > 0 (neglecting  the  trivial  case  ^ = 0)  by 
construction  of  the  algorithm  and  the  definition  of  ^ • Now  using  the 
equation  in  step  3,  one  finds  that 

CTK+ 1 ctk  " “ l "NA'Cl  + € ) ) (°K  ' ^1 } 

(4. 4.11 


where  1 > r = a / (N"A(1  + 0)  > 0.  This  implies  that  a converges  to 
at  least  geometrically,  and  the  convergence  of  AL-(u,  £,y)  is  proven. 
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A lower  bound  on  $ J Let  S be  any  set  of  arcs  in  the  network. 

S should  be  chosen  as  an  estimate  of  the  bottleneck  associated  with  the 
linear  routing  problem  to  achieve  a good  lower  bound  on  p^.  Further- 
more, define: 


P S(CT,  f) 


i<§ 


1 - p.'a 

i 


Using  the  same  methodology  employed  in  Section  2.1,  one  can  find  a 
lower  bound  for  the  value  of  Pg(a,  g)  for  any  geF.  Therefore, 

€g(a,f)  a <7Ps(ff,f),  f - <p)  * Ps(o,f)-Ps(ff,g) 

where  <p  solves  min  (VP Q(cr,  f),  h).  In  particular,  let  g be  any  solution 
heF 

to  4. 1.  1.  Then  it  is  clear  that 


Ps«Ug>  £ |S|  — x 

1 - \ /a 


since  every  component  of  g satisfies  f./Ch  s p^.  Thus, 


Is  | I * Ps(cr,g)  * Pg(cr,f)  - Co  (a,  f) 

1 - ^ /a  b s s 


which  yields  a lower  bound  for  $ that  is;  , 


^1  * a Ps(a,  f)  - Cg(o,  f)  = ^(S.a.f) 


Some  trade-offs.  Now  suppose  that  f^  approximately  solves 

min  P(a  g)  for  a aK  "very  close"  to  and  an  c sufficiently  small. 

gCF  1 

Then  S should  be  easy  to  identify,  since  channels  which  can  have  a 


t For  a proof  of  the  convergence  of  the  lower  bound,  see  Appendix  A, 

Thm  4.  4.  A. 


utilization  factor  ”farM  from  or  will  tend  to  be  driven  down  by  the 

K. 

objective  function.  Also,  eg(c_,  f ) should  be  small  since  the  terms 
of  Pc(al_,  f)  are  the  dominating  terms  of  P (a, . , f)  in  the  neighborhood  of 

o K.  tv 

K K ^ 

f . Under  these  assumptions  ft  (S,cr. f ) will  not  be  far  from-ft-<as 

one  can  see  by  the  substitutions  e0(cr  f*x)s=  0 and  Pc(a  f^ ) =3 
- b h.  b k 

I S 1/  (1  into  the  above  expression  for  /y(S,oy.,  f^)). 

Some  trade-offs  are  now  presented.  One  cannot  make  e 
arbitrarily  small  because  at  some  point  the  rate  of  convergence  of  the 
AL*(a ,C,y)  algorithm  will  suffer,  and  yet  a fairly  small  ( is  clearly 
necessary  for  the  lower  bound  on  ^ to  be  tight. 

A "lower  bound  directed"  algorithm  would  keep  e small.  One 
might  embellish  the  algorithm  by  heavily  weighting  those  terms  of 
P(ct  g)  corresponding  to  the  channels  which  one  estimates  are  most 
likely  in  the  bottleneck.  If  one  guesses  right,  the  behavior  of  P(ft^,  f) 
in  the  neighborhood  of  f^  will  be  more  similar  to  that  of  Pg(a^.,  f),  and 
the  control  on  the  magnitude  of  es(cr  f ) will  be  enhanced. 

Despite  the  foregoing  possibility,  a "convergence  directed" 
algorithm  was  found  to  be  more  useful  for  the  purposes  of  this  thesis 
after  experimentation  with  each  type.  In  essence  one  opts  to  allow 
the  lower  bound  to  degrade  in  an  attempt  to  speed  up  the  algorithm. 
Towards  this  same  goal,  minimization  of  P(a,  f)  was  replaced  in  the 
definition  of  AL* (a , € , y)  by  a minimization  of 

P*(a,f)  = 2 m 'vFC-— ) 

1-1  11 


The  advantage  of  P*(a,  f)  over  P(cr,  f)  lies  in  the  form  of  the  gradient 


-oo- 


k 


whose  terms  (i.e.  OP- (<J,  f)/£>f.  = l/(aC.  -f  ))  are  first  rather  than 

second  order  poles.  Due  to  the  number  of  times  which  the  gradient 

is  calculated,  the  simplicity  of  the  gradient  translates  into  significant 

time  reductions.  The  disadvantage  is  that  the  behavior  of  P *(o  f)  in 

the  neighborhood  of  f1'"  does  not  resemble  that  of  PQ(<V,  0 as  well  as 

the  behavior  of  P(^K>  f)  does.  This  tends  to  degrade  ^(a^,  f1'")  and 

consequently  0 (S,cr„,  f^). 

Li  K. 

The  deterioration  of  the  lower  bound  is  in  fact  so  severe  that  its 

utility  is  virtually  destroyed.  One  is  therefore  motivated  to  examine 

how  0^(S,(J,  f^)  behaves  as  a function  of  O ; specifically,  no  longer  is 

the  cr  of  0T  (S.cr,  f ) taken  to  be  synonymous  with  the  cr  of  P*(ct  f). 
l*  K.  K 

As  CT  grows  larger  than  the  behavior  of  Pg(cr,  f)  in  the  neigh- 

borhood  of  f becomes  even  less  similar  to  that  of  P*(ff  , f),  which  is 

K. 

discouraging.  On  the  other  hand,  as  a grows  and  Pg(a,  f)  tends  to 
behave  more  linearly  about  f^,  cg(a,  f1^)  becomes  a less  liberal 
estimator  of  the  maximal  error  in  Pg(cr,  f),  which  favors  the  growth  of 
cr.  Consider  the  following: 
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= lim 
a ■*  *> 


(CT 


I a 

igS  ’N 

l-2p./j  4-  / 

LI  11 

PJ 

<PcF  ifS  (i  -p  /a) 


lim 

a -*  ® 


, a - 2 p + o.  / C . v 

v*  i i 1 ' _ v ,-t 

min  .If  9 ; • -t  u 

<PeF  ucS  (1  . p /a)“  ' i«S 

1-2  p. /a  + (Pj/ac. 

min  E ^ 

<Pc F i«S  (1  . p/cr)^ 


min 


<Pj/Ci  * P j/a 


<PcFicS  (i-p/a)" 


lim 

a -*  cd 


min 


1 -2  p. /a  + a/ac . 


pfF  US  (i  -p./a)^ 


min  E0  (<p. /C.) 
<PcF  ifb  i 1 

|s| 


In  order  to  gain  insight  into  the  utility  of  /^(S,®,  f),  the  proLm.  m 
can  be  simplified.  Let  C.  = 1 for  every  ie{l,  2,  • • • , NA] , so  that 

^(S,-,f)  . (1/lsDmin.Egtfj 

This  can  be  evaluated  via  a shortest  route  algorithm  with  a distance  of 
1 assigned  to  each  arc  which  is  estimated  to  be  in  the  bottleneck  and  a 
distance  of  0 on  the  remaining  channels. 

A 

The  condition  for  which  ^(S,  f)  is  in  fact  ^ is  derived  as  follows. 
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Let  g(F  be  a solution  of  4.  1.  1 , and  assume  that  g is  known  as  a convex 


N N mk 


- • rvi  Lr 

combination  of  chains,  i.e.  g E E E X ■ , . For  S 

k = l m = l j=l  Jmk  (j) 
m/k 

chosen  as  the  bottleneck,  it  follows  that; 


k = (1/  IS  | ) E g 

1 1 £ S 1 


mk 


. . N N cmk 
(1/  S ) E E E , r t i~ 

k = 1 m = l j = l Ajmk  mk  (j) 
mrk 


where  = the  number  of  bottleneck  channels  used  by  chain 
Thus; 


(1/  |S|)  S L r i*J  flkX.  ' N 

k=l  m=  1 mk  mk\j  = i Ajmk/ 

m^k 


+ ''dsl'l  L jCk  xjmk  (X—  - 

m^k  J 


where  j^^is  the  minimum  number  of  bottleneck  channels  on  any  path 
from  m to  k.  This  yields: 

K = f)  + Dif(g) 

where  Dif(g)  is  defined  to  be  the  latter  term  in  the  equality  previous  to 
its  appearance.  Thus  Dif(g)  = ®,  f)  and  is  independent  of  g 

(providing  g solves  4.  1.  1).  This  yields  the  desired  condition; 


“ GO  “ 


If  any  gcF  which  solves  4.  1.  1 is  such  that  every  active 
associated  chain  uses  a minimal  number  of  channels 
which  are  members  of  the  bottleneck,  then  the  lower 
bound,  j3^(S,  f),  will  be  equal  to  when  the  choice  of 

S is  correct. 

Consult  Figs.  4.3.1  and  4.3.2,  both  of  which  display  very  simple 
cutsets  as  the  bottleneck.  It  is  clear  that  f)  can  be  made  equal 

to  for  both  of  these  networks  since  any  commodity  which  must  use 
the  bottleneck  will  pass  through  no  more  than  one  bottleneck  channel, 
and  all  other  commodities  will  avoid  the  bottleneck  entirely. 

Figure  4.3.  3 displays  a situation  where  it  is  not  clear  whether 
the  above  condition  holds.  However,  in  the  solution  of  the  associated 
linear  routing  problem  generated  by  the  LCH  algorithm,  the  following 
occurs; 

(i)  The  minimum  number  of  bottleneck  channels  between 
nodes  4 and  7 is  one,  e.g.  as  evidenced  by  path 
(4,  5,  7). 

(ii)  The  chains  used  for  the  commodity  (4,  7)  in  the 
solution  of  4.  1.  1 generated  by  the  LCH  algorithm 
were: 

92.  6%  on  path  (4,  5,  7 ) 

1. 6%  on  path  (4,  6,  7) 

5.  8%  on  path  (4,  2,  7 ) 

(iii)  Both  of  the  latter  paths  are  active,  and  they  use 
more  than  the  minimal  number  of  bottleneck 


f 


, 


channels  between  nodes  4 and  7. 


Thus  f)  cannot  be  made  as  great  as  p when  S is  chosen  as  the 

bottleneck,  and  the  required  difference  could  be  evaluated  if  desired 
according  to  the  definition  of  Dif(g)„ 

It  follows  that  how  /3^(S,a,  f)  varies  as  a function  of  O’,  and  in 
particular  which  value  of  o provides  the  largest  lower  bound,  depends 
on  the  problem  at  hand.  Many  times,  regardless  of  whether  the 
AL*(Q! , €,  y)  algorithm  is  "lower  bound  directed"  or  "convergence 
directed",  the  best  value  of  o occurs  at  a = °°  where  the  bound  is 
exact  for  a properly  chosen  S.  Often  ^(S,00,  f)  (for  S chosen  as  the 
bottleneck)  is  not  exact  but  not  far  below  j^,  in  which  case  the  "lower 
bound  directed"  algorithm  might  find  that  a o near  a provides  the 
best  bound,  while  the  "convergence  directed"  algorithm  is  more 
likely  to  find  that  a - <*>  is  best.  More  will  be  said  on  this  subject  in 
the  next  section. 

4.5  Computational  Experience  with  the  AL';:(tt,  g,y)  Algorithm 

The  computational  results  of  this  section  are  important  for 
more  than  merely  the  AL*(a,  e,y)  algorithm.  This  section  will 
confirm  the  speculations  of  Chapter  II  on  the  relative  performance  of 
the  EF,  TR,  and  CH  algorithms.  To  this  end,  the  same  four 
examples  which  were  under  consideration  in  Section  4.3  are  considered 
again  in  this  section. 

The  selection  of  parameters.  The  parameters  a,  c,y  of  the 
AL*(a ,<,y)  algorithm  are  relatively  unimportant  unless  one  makes 
immoderate  choices  for  them.  Experience  with  the  algorithm 
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indicates  that  the  successive  problems  solved  in  step  2 of  the 
algorithm  should  not  be  drastically  different,  and  yet  should  d;  ffe: 
enough  so  that  the  step  has  some  utility.  This  suggests  that  > she 
be  neither  particularly  close  to  its  lower  bound  of- 1 nor  extremeh 
large;  using  the  same  reasoning,  a should  not  be  near  its  lower 
bound  of  0 nor  its  upper  bound  of  1.  Also  e should  not  be  chosen 
unreasonably  small,  else  the  minimization  of  step  2 will  be  needle  o 
accurate.  A detailed  description  of  the  experimentation  with  these 
parameters  is  consequently  of  little  importance.  It  suffices  to  sa> 
that  experimentation  led  to  the  chosen  values  of  > = 5,  a = 0.  6,  and 
€ = 0.  16.  From  this  point  forward,  the  parameter  list  will  be 
dropped  from  the  names  of  the  algorithms,  with  the  implication  that 
the  above  values  are  employed. 

Relationship  to  the  TR,  EF,  and  CH  algorithms.  Later  in  tin  - 
section,  the  TR*,  EF  , and  CH*  algorithms  will  be  compared  almo  t 
as  if  they  were  the  TR,  EF,  and  CH  algorithms.  There  are  several 
reasons  and  justifications  for  this  comparison.  Since  the  TR,  Eh, 
and  CH  algorithms  all  solve  the  minimum  delay  problem,  there  is  tin- 
difficult  question  of  how  heavily  the  network  should  be  loaded  (or  how 
the  requirements  should  be  scaled)  for  a meaningful  comparison. 
Obviously,  if  the  demand  on  the  network  is  extraordinarily  low,  the 
problem  can  virtually  be  solved  with  one  shortest  route  computation. 
This  extreme  is  not  desirable  for  comparison  purposes;  however,  if 
the  demand  is  as  great  as  the  network  will  allow,  then  the  algorithm  ■ 
will  be  compared  near  their  "worst",  with  nearly  as  many  columns  .n 
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the  problem  as  any  routing  problem  will  ever  require  for  that  network. 

The  latter  extreme  is  precisely  that  which  the  TR  . EF*.  and  CH* 
algorithms  seek.  Furthermore,  the  form  of  the  objective  function 
used  in  the  TR*.  EF*,  and  CH*  algorithms  is  similar  in  character  to 
that  of  the  expected  delay  function  of  Chapter  II;  i.o.  each  goes  to 
infinity  as  the  flow  in  any  channel  approaches  the  capacity  of  the 
channel.  In  essence,  then,  the  EF*,  TR",  and  CH  : algorithms 
solve  a sequence  of  problems  very  similar  to  minimum  delay  problems 
over  a variety  of  loads  on  the  network. 

It  is  also  noted  that  using  the  EF*  algorithm  is  an  extremely 
attractive  way  of  starting  the  EF  algorithm  (and  similarly  for  the 
CH  and  TR  algorithms).  At  least  it  is  the  best  of  the  methods  that 
this  experimenter  has  attempted.  When  is  driven  below  unity  in 
the  EF*  algorithm,  then  the  EF  algorithm  has  a feasible  (and  reason- 
able) starting  point.  When  the  EF  algorithm's  objective  function  is 
complicated,  unexpected  benefits  may  occur.  A great  deal  of  effort 
may  have  been  saved  by  using  the  objective  function  of  the  EF* 
algorithm  until  the  estimate  of  the  solution  is  in  the  correct  ballpark. 

Discussion  of  the  examples.  The  examples  are  presented  in 
turn,  along  with  a discussion  of  each.  With  each  set  of  curves,  the 
curve  of  the  fastest  linear  program  from  those  in  Section  5.  3 is 
included.  Every  algorithm  had  the  same  initial  flow,  although  that 
may  not  be  evident  since  the  curves  are  not  shown  until  they  are  in  the 
range  of  the  graph. 

I 


Figure  4.  5.  1 displays  the  curves  for  the  linear  routing  problem 


Fig.  4.  5.1  The  experimental  convergence  of  /S ^ of  the  linear  routing 
problem  associated  with  Fig.  4.  3.  1 for  the  LTR,  TR*. 
and  FW algorithms.  If  the  curves  for  the  CH*  and  EF': 
algorithms  had  been  drawn,  they  would  nearly  coincide 
with  that  of  the  TR  algorithm.  The  curves  under  the 
dotted  line  are  the  associated  lower  bounds  on  as  a 
function  of  time. 
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associated  with  the  ten  node  network  of  Fig.  4.  3.  1.  This  networx  has 
30  directed  arcs  or  an  average  of  3 outgoing  arcs  per  node  and  33%  of 
the  number  of  arcs  of  a fully  connected  network.  The  FW*  algorithm 
was  much  faster  for  this  example  than  the  EF'  , TR*.  and  CH* 
algorithms,  although  this  was  not  typical  of  the  other  examples.  The 
following  is  a list  of  the  number  of  columns  and  the  number  of  storage 
locations  required  at  the  final  estimate  of  the  solution  (including 
routing  matrices  for  the  EF*  algorithm): 


# of  columns 

# of  locations 

EF* 

17 

2210 

TR* 

45 

900 

CH* 

142 

369 

It  is  noted  that  the  TR*  and  CH*  algorithms  exceed  their  "bounds"  of 
Chapter  II  on  the  maximum  number  of  columns,  although  one  could 
pivot  out  the  excess  columns  if  desired. 

The  lower  bound  on  /3  became  exact  at  some  point  in  each  of  the 
EF*,  CH*.  TR*,  and  FW*  algorithms  for  the  10  node  problem.  This 
is  not  to  imply,  however,  that  the  greatest  lower  bound  could  not  have 
been  achieved  earlier.  Much  depends  on  the  method  of  choosing  S, 
and  one  might  try  several  different  methods  of  choosing  S each  time 
the  lower  bound  is  calculated.  However,  a great  deal  of  experimenta- 
tion with  methods  for  choosing  S has  not  been  done.  For  these 
examples,  only  one  method  was  used,  and  it  had  the  property  of  being 
rather  conservative  (i.e.  tending  to  include  too  many  channels  in  S) 
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early  in  the  run,  and  less  so  as  the  algorithm  progressed.  The 
intention,  which  normally  was  realized,  was  to  obtain  a good  lower 
bound  at  some  point  in  the  run  of  the  algorithm  without  spending  a 
great  deal  of  computation  time  on  it.  For  these  network  examples, 
all  of  which  were  solved  via  "convergence  directed"  algorithms,  the 
best  lower  bound  on  $ was  invariably  found  at  o = *>. 

Figure  4.5.2  displays  the  curves  for  the  linear  routing  problem 
associated  with  the  26  node  network  of  Fig.  4.3.2.  This  network  has 
60  directed  arcs  or  an  average  of  2.  31  outgoing  arcs  per  node  and  9% 
of  the  maximum  number.  The  remarks  of  Chapter  II  on  the 
comparisons  of  the  algorithms  indicate  that  the  CH  ■ algorithm  would 
not  be  as  effective  as  the  other  algorithms  on  a very  sparsely 
connected  network,  and  this  is  reflected  in  the  curves. 

The  error  bound  was  the  least  effective  on  this  problem,  although 
the  bottleneck  is  simple  (see  Appendix  E).  From  Fig.  4.3.2,  one  can 
see  that  arcs  (14,  15),  (20,25),  and  (24,4)  form  the  bottleneck.  It 
also  happens  in  the  solution  that  the  utilization  factor  of  arc  (15,  16) 
turns  out  to  be  nearly  as  great  as  and  the  programs  normally 
chose  all  four  arcs  to  be  the  estimate  of  the  bottleneck.  With  so  few 
arcs  actually  in  the  bottleneck,  the  resultant  percentage  error  in  the 
lower  bound  from  the  true  value  of  $ is  large.  A bit  of  good  fortune 
went  into  play  for  the  FW  algorithm,  as  it  once  (incorrectly)  chose 
arcs  (15,  16),  (20,25),  and  (24,4)  as  the  bottleneck,  which  nevertheless 
yielded  a fairly  good  lower  bound  on  $ . When  there  are  a larger 
number  of  arcs  in  the  bottleneck,  ^(S,a,  g)  is  not  so  sensitive  to  the 
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Fig.  4.  5.  2 The  experimental  convergence  of  ^ of  the  linear  routing 
problem  associated  with  Fig.  4.  3.  2 for  the  LTR,  EF*. 
CH^,  and  FW*  algorithms.  If  the  curves  for  the  TR* 
algorithm  had  been  drawn,  they  would  nearly  coincide  with 
that  of  the  EF*  algorithm.  The  curves  under  the  dotted 
lines  are  the  associated  lower  bounds  on  % as  a function 
of  time. 
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choice  of  S. 

Figure  4.  5.  3 corresponds  to  the  problem  associated  with  the 
seven  node  symmetric  network  depicted  in  Fig.  4.  3.  3.  Note  that  for 
the  first  time  the  "knee"  of  every  convergence  curve  of  each  AL. 
algorithm  occurs  before  the  knee  of  the  linear  programming 
algorithm.  This  network  is  the  most  fully  connected  under  considera- 
tion, with  66%  of  the  number  of  channels  of  a fully  connected  network. 

It  is  therefore  noted  without  surprise  that  the  performance  of  the  CfT 
algorithm  timewise  was  better  than  the  TR*  and  EF*  algorithms. 

The  CH*  algorithm  required  only  61  chains  (most  of  which  have  only 
one  or  two  elements),  the  TR*  algorithm  required  31  trees  (each 
having  6 elements)  and  the  EF*  algorithm  used  21  extremal  flows  (each 
having  28  elements  plus  a routing  matrix  of  49  elements).  This  was 
the  first  network  such  that  ^ (S,  ®,f)  could  not  equal  $ for  any  choice 
of  S,  but  the  lower  bound  differed  from  ^ by  less  than  £%. 

Finally,  consider  Fig.  4.5.4  which  corresponds  to  the  problem 
associated  with  the  eleven  node  symmetric  network  shown  in  Fig.  4.  3.  4. 
The  knee  of  the  convergence  curves  of  the  AL*  algorithm  is  even 
further  to  the  left  of  the  knee  of  the  linear  programming  curve.  This 
network  is  less  fully  connected  than  the  previous  one,  having  40%  of 
the  channels  of  a fully  connected  network,  and  the  EF* , CH*.  and  TR 
algorithms  performed  nearly  equally  timewise.  At  the  final  estimate 
of  the  solution,  the  CH*  algorithm  had  161  chains  (most  of  which  had 
no  more  than  3 or  4 elements),  the  TR:  algorithm  had  59  trees  (each 
having  10  elements),  and  the  EF  algorithm  had  39  extremal  flows 
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The  experimental  convergence  of  Si  of  the  linear  routing 
problem  associated  with  Fig.  4.3.3  for  the  LCII,  EF;:, 
CH*,  and  FW * algorithms.  If  the  curves  for  the  TR  ! 
algorithm  had  been  drawn,  they  would  be  in  between  those 
of  the  EF*  and  CH*  algorithms.  The  curves  under  the 
dotted  line  are  the  associated  lower  bounds  on  $.  as  a 
function  of  time. 


Fig.  4.  5.  4 The  experimental  convergence  of  3^  of  the  linear  routing 
problem  associated  with  Fig.  4.3.4  for  the  LCH*,  CH';:, 
and  FW  algorithms.  If  the  curves  for  the  TR*  and  EF* 
algorithms  had  been  drawn,  they  would  nearly  coincide 
with  that  of  the  CH;:  algorithm.  The  curves  under  the 
dotted  line  are  the  associated  lower  bounds  on  3 
function  of  time. 
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(each  having  -14  elements  plus  a routing  matrix  of  121  elements). 

Once  again,  the  TR*  and  CH"  algorithms  exceeded  their  "bounds''  on 
the  bumber  of  columns,  but  not  by  much.  Although  j^(S,  <*>,  f)  could 

A 

not  be  made  exact,  lower  bounds  about  2%  beneath  the  value  of  8 ^ were 
achieved. 

As  the  problems  became  more  complicated  (in  terms  of  the 
number  of  binding  constraints  in  the  solution),  the  trend  is  very  clear. 
The  "knee"  of  the  convergence  curves  of  the  AL*  algorithms  occurs 
far  before  the  linear  program  converges  to  the  solution.  If  one  is 
content  with  the  inaccuracy  that  exists  at  the  knee  of  the  convergence 
curve  (which  is  less  than  about  2%  for  these  problems)  and  one 
expects  the  problem  to  be  complicated,  then  the  nonlinear  approach 
may  be  the  preferred  approach.  If  one  wanted  to  run  the  programs 
only  until  the  knee  of  the  convergence  is  reached,  then  it  is  probably 
in  his  best  interest  to  try  several  means  of  choosing  S at  the  termina- 
tion in  order  to  get  the  best  lower  bound  on  the  objective. 

One  practical  advantage  of  the  nonlinear  approach  is  that  it  also 
tends  to  be  searching  for  a solution  which  takes  into  account  higher 
order  criteria  for  the  solution.  For  example,  in  the  LTR  solution  to 
the  problem  associated  with  Fig.  4.  3.  2,  there  are  ten  active 
constraints,  although  oniy  three  are  required  to  be  active.  As  one 
can  see  from  the  form  of  the  objective  functions  for  the  nonlinear 
approach,  those  utilization  factors  which  are  not  necessarily  at  ^ will 
tend  to  be  driven  to  a lower  value.  Thus  although  the  estimate  of  the 
solution  to  the  linear  routing  problem  which  the  nonlinear  approach 
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provides  for  the  first  criterion  has  a higher  value  for  the  objective 
than  the  linear  program,  other  aspects  of  the  solution  as  viewed  from 
the  higher  order  criteria  may  be  better  (see  Appendix  D). 


CHAPTER  V 


STATISTICAL  LOAD  SHARING  IN  A COMPUTER- 
COMMUNICATION  NETWORK 

This  chapter  deals  with  the  application  of  the  concepts  developed 
in  this  thesis  to  the  problem  of  load  sharing  in  a computer-communication 
network.  A set  of  computers  is  assumed  to  be  connected  via  a store 
and  forward  network.  Each  computer  receives  jobs  which  it  must 
process  and  return  to  the  customer.  Since  the  computers  are  inter- 
connected, a heavily  loaded  computer  has  the  option  of  sending  jobs 
elsewhere  via  the  communication  network,  having  the  jobs  processed  at 
that  location,  and  having  the  results  returned  via  the  network.  In  this 
way,  the  resources  of  the  entire  system  can  be  pooled  for  the  overall 
benefit.  The  model  of  this  system  which  is  employed  in  this  chapter 
is  that  of  Wunderlich  [WUND75]  and  McGregor  [McGR74]. 

This  chapter  is  not  directed  towards  providing  a large  amount  of 
qualitative  insight  on  the  load  sharing  problem.  The  general  subject 
matter  of  load  sharing  is  handled  in  greater  depth  in  Wunderlich' s work 
[WUND76];  however,  the  details  of  the  algorithm  which  is  applied  to 
the  statistical  load  sharing  problem  are  more  complete  in  this  thesis. 

5.1  Models  and  Definitions 

The  model  of  a computer  used  in  this  study  is  an  M/M/1  queue. 

This  is  a simple  model  of  a computer  operating  in  a batch  processing 
mode  [KLEI75],  For  this  model  it  is  assumed  that: 
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1.  The  input  stream  of  jobs  to  the  jth  computer  from 
outside  the  system  is  a Poisson  process  with  mean 
arrival  rate  Xj  jobs  per  unit  time. 

2.  The  number  of  operations  required  per  job  is  dis- 
tributed as  a negative  exponential  with  mean  III. 


The  jth  computer  performs  FT  operations  per  unit 
time.  This  means  that  the  service  time  per  job 
(not  including  waiting  time)  is  distributed  as  a 
negative  exponential  with  mean  1/iR  , 


4.  Jobs  are  processed  in  a first-come-first-served 

manner.  If  the  computer  is  busy  when  a job  arrives, 
it  is  queued  in  an  infinite  buffer. 


The  model  used  for  a store  and  forward  communication  channel 
is  also  an  M/M/1  queue.  For  the  communication  channel  it  is 
assumed  that: 


1.  The  length  of  messages  in  bits  is  distributed  as  a 
negative  exponential  with  mean  1 for  programs 
(computer  inputs)  and  mean  1/Mr  for  results 
(computer  outputs). 

2.  The  ith  communication  channel  has  a capacity  C. 
bits  per  unit  time  so  that  the  time  to  transmit  a 
message  is  also  distributed  as  a negative  exponential 
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with  mean  l/H  C.  or  1 //i  (J.. 

'p  x r 1 

3.  If  a message  passes  through  more  than  one  com- 
munication channel,  its  length  is  assumed  to  be 
chosen  independently  at  each  channel  through  which 
it  passes. 


This  queueing  model  of  a store  and  forward  communication  channel 
has  been  extensively  studied  by  Kleinrock  [KLEI64]  . 

One  final  assumption  that  must  be  made  is  that  the  message 
lengths  of  programs,  the  number  of  computer  operations  they  require 
and  the  message  lengths  of  the  results  are  all  independent  random 
variables.  This  assumption  is  required  to  make  the  problem 
mathematically  manageable. 

The  performance  measure  used  in  this  study  is  the  steady  state 
expected  time  to  process  a computer  job.  For  a system  of  N 
computers,  this  expected  time  is  given  by: 


T 


X. 


j = 1 \r 


T. 

J 


(5.  1.1) 


N 

where  = L 

1 3 = 1 J 


The  time  Tj  is  the  expected  time  to  process  a computer  job  that  enters 
the  system  at  the  jth  computer.  This  time  includes  the  computation 
time  required  by  the  job  and,  if  the  job  is  processed  by  a computer 
other  than  the  one  at  which  it  was  submitted,  it  also  includes  the  com- 
munication  time  required  to  send  the  program  to  the  processing 
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computer  and  the  time  to  return  the  results  to  the  point  of  origin. 

5.  2 Formulation  of  the  Load  Sharing  Problem  as 
a Multicommodity  Flow  Problem 

This  section  considers  the  formulation  and  solution  of  the  statis- 
tical load  sharing  problem  as  a multicommodity  flow  problem. 
Presentation  of  the  formulation  is  best  done  via  a small  example.  A 
three-computer  network  is  depicted  in  Fig.  5.2.  1.  The  network 
topology  for  representation  as  a multicommodity  flow  problem  is  shown 
in  Fig.  5.2.2.  where  each  computer  is  associated  with  two  nodes  (an 
input  and  an  output  node),  and  each  physical  communication  channel  is 
represented  by  one  channel  for  programs  and  a second  for' results. 
These  following  definitions  are  required  for  the  multicommodity  formu- 
lation: 

fj  = average  demand  on  computer  j (jobs  per 
unit  time). 

= average  demand  on  channel  i by  jobs  in 
program  form  (jobs  per  unit  time). 

fri  = average  demand  on  channel  i by  jobs  in 
result  form  (jobs  per  unit  time). 

NA  = number  of  communication  channels. 

Since  the  only  delay  associated  with  the  communication  channels 
is  assumed  to  be  the  queueing  delay,  the  average  delay  incurred  by  a 
job  (in  program  or  result  form)  on  channel  i may  be  approximated  bv; 
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Fig.  5.2.1  A three  computer  network  connected  by  a 
communication  network. 


Fig.  5.2.2  The  representation  of  the  system  of  Fig.  5.2.  1 
for  a multicommodity  flow  problem. 
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which  agrees  with  the  vv on  -k 
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■age  delay  formula  for  the  .U/'M/l 


queue  model  of  a communication  channel  in  the  case  f . = 0 or  f . = 0 
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or  fi  = ' The  average  delay  of  a randomly  selected  job  passing 


through  a system  of  computers  is  then  approximated  by: 
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which  is  a consequence  of  5.  1.  1,  5.2.  1.  and  the  fact  that  the  entire 


system  has  been  modeled  as  a network  of  queues  f.JACK57j. 


Define  f . to  be  the  average  demand  on  computer  j due  to  jobs 


which  enter  and  leave  the  system  at  computer  k.  That  is. 


N v k if 

f.  = £ f.  . In  an  analogous  manner,  also  define  f . and  t ..  Then 
J k = 1 J ° pi  ri 


for  every  j,  k € Cl , 2,  • • • , N}  , one  may  write  two  conservation  of  flow 
equations  (which  are  most  easily  understood  in  conjunction  with 
Fig.  5.2.2); 


Discussion  of  this  approximation,  and  of  additional  conditions  under 
which  Equation  5.2.1  is  exact,  is  deferred  to  Section  5.5. 
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where  the  computers  of  the  origin  and  destination 

of  the  ith  channel 

are  OR(i)  and  DS(i)  respectively.  Let  fk  = [fk 

^ Lpl  p2  ’ pNA  rl* 

^2'  * * * * ^NA’  fi*  * * * * fjsr]  » and  define  f = jS  The  vector  f is 

said  to  be  an  element  of  the  set  F if  there  exists  a collection, 

such  that  (i)  f^  ^ 0 for  every  k,  (ii)  f^  satisfies  the 

above  conservation  of  flow  equations  for  every  j and  k,  and 
N Jc 

(iii)  £ i = f.  Then  the  routing  problem  may  be  stated  as; 
k = l 

minimize  T (f)  (5.2.3) 

feF 

where  T(f)  A « if  (f  ,/ji  + f ./ji  ) ^ c.  for  anv  i or  f.  >£R.  for  anv  i. 

= pi  'p  ri  'r  i 7 j j ~ 

U HU 

If  f solves  5.2.3  and  r = 0 or  f = 0 for  every  ie{l,  2,  • • • , NA] , then 

pi  ri  J 9 

u 

the  expression  for  the  average  delay  becomes  exact  at  f subject  to  the 
conditions  of  the  model.  For  many  practical  problems  this  will  be 
the  case;  however,  the  many  problems  involving  links  which  are  not 
duplex  and  some  problems  with  all  duplex  links  (see  Fig.  5.2.3),  the 
expression  for  delay  remains  an  approximation. 
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5.3  An  Efficient  Algorithm  to  Solw  the  General 

Computer  Load  Sharing  Problem 

At  this  point,  an  objective  function  has  been  derived  which  is 
convex  and  has  a gradient  with  all  components  non-negative  over  F so 
that  a variety  of  algorithms  may  be  directly  applied  to  solve  5.2.3. 
Decomposition  techniques,  such  as  the  extremal  flow  (EF)  algorithm 
of  Cantor  and  Gerla  [CANT7  4]  are  particularly  attractive.  This 
section  presents  modifications  to  the  EF  algorithm  which  produce  an 
efficient  algorithm  for  the  load  sharing  problem.  Initially,  it  is  noted 
that  5T(f)/0fr.  = (^p/jir)(0T(f)/&fp.),  and  this  fact  is  exploited  in  both  the 
master  problem  and  the  subproblem  of  the  EF  algorithm  to  reduce  the 
number  of  arcs  and  the  number  of  nodes  which  are  considered 
explicitly  by  the  algorithm. 

The  master  problem  of  the  EF  algorithm  is: 

~ / M /;  \\ 
minimize  T [ L a.c p ) 

^i=  1 1 / 

M 

subject  to  L a.  = 1,  and  a.  ^ 0 
i=l  1 i 


•where  \ is  a generated  set  of  extremal  flows  of  the 

set  F.  The  fundamental  operation  necessary  to  carry  out  the  mini- 


mization is  the  calculation  of  — — T 1 

0 a. 

J 


fL  a<P(i>V 

' i-1  i J 


If  VT(f)  = 


!?T(f)/ofpl.oT(f)/8fp2,  • • • ,af(f)/ofpNA,aff)/afrl.of(f)/6fr2,  • • • . 

vT(f)/5fr^^,jjT(f)/0fj  .aTlfj/Ofg,  * • * , STIO/Ofj^)  , (,  ) denotes  the  usual 

inner  product  in  Euclidian  space,  and  f*  = 2?  asf>^  (where  f is  the 

i=1  1 


current  estimate  of  f ),  then 
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computation  of  the  subproblem  may  be  done  first  over  the  networn  of 
primed  nodes,  then  over  the  unprimed  nodes,  and  finally  combined. 
Further,  the  network  distances  between  the  primed  nodes  only  differ 
by  a multiplicative  constant  from  the  network  distances  between  the 
corresponding  unprimed  nodes,  so  that  the  shortest  routes  may  be 
found  over  the  unprimed  network  with  the  shortest  routes  over  the 
primed  network  being  implicit.  Ultimately,  only  the  shortest  routes 
between  j and  j1  for  every  jc { 1 , 2,  • • • , Nj  are  required  (since  they  are 
the  only  origin-destination  pairs  with  a nonzero  demand  for  a flow  of 
jobs),  and  each  may  be  found  via  the  solution  of 


minimize  <d.k  + dkk,  + d^.,: 


(5.  3.  1) 


1 s k s N 


JU  « ( f 4 

where  d^  is  the  distance  of  the  shortest  j to  k path,  d^,  — - 


jk 


and  d'  = (ji  j a )d,  ..  Thus  the  number  of  nodes  explicitly 
k 3 p 'r  tvj  r 

represented  is  reduced  by  a factor  of  two  for  the  subproblem.  As 
shortest  route  computation  time  may  grow  as  fast  as  the  number  of 


nodes  cubed,  this  represents  up  to  a factor  of  eight  reduction  in 

vr2 


shortest  route  computation  time  (less  the  order  N operations  of 
5.  3.  1,  which  generally  should  be  small  in  comparison  to  the 
reduction). 

The  final  modification  to  the  EF  algorithm  is  to  base  the 
decomposition  on  chains  rather  than  extremal  flows  (each  extremal 
flow  is  a superposition  of  N chains).  This  provides  greater  freedom 
in  finding  downhill  directions  resulting  in  fewer  shortest  route 


r 


L 
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computations,  and  memory  requirements  are  reduced  as  chains  can 
be  stored  compactly.  The  details  of  this  alteration  were  indicated  in 
Chapter  II.  This  load  sharing  algorithm  based  on  a chain  decom- 
position is  denoted  the  LSCH  algorithm. 

In  conclusion  of  this  section,  the  following  remarks  contrast  the 
LSCH  algorithm  with  that  of  McGregor  and  Boorstyn  [McGR75]  who 
devised  the  only  other  known  algorithm  for  the  problem: 

(i)  The  LSCH  algorithm  is  free  of  the  constraint  that 
results  must  travel  the  reverse  route  of  the 
respective  program.  This  leads  to  lower  values 
of  the  objective  and  allows  a more  general 
communication  network. 

(ii)  An  upper  bound  to  the  absolute  error  of  any 

estimate  of  the  solution  to  5.2.3  is  a by-product 
of  the  LSCH  algorithm,  which  allows  control  over 
the  necessary  running  time  of  the  algorithm  for 
any  desired  accuracy. 

(iii)  The  routing  strategy  is  not  derived  from  the 
solution,  but  is  explicitly  the  solution  itself, 
i.e.  , the  weight  of  each  chain  is  the  fraction  of 
jobs  arriving  at  the  source  of  the  chain  which 
should  take  the  route  described  by  the  chain 
through  the  network. 


L . . A 


-n  3- 


5 . 4 An  Example  of  Load  Sharing 

A small  example  of  the  use  of  this  algorithm  is  presented.  The 
example  is  of  a ten-computer  system  for  each  of  which  R.  = 10  steps / 

7 

sec.  and  l - 4 x 10  steps/ job.  The  communication  network  is  shown 

in  Fig.  5.4.1  and  consists  of  30  directed  channels,  each  having 

capacity  9.6  Kb /sec.  The  mean  lengths  of  messages  are  1/M  =50  Kb 

P 

and  = 500  Kb.  A randomly  generated  set  of  input  rates  to  the 
computers  with  and  without  load  sharing  are  given  in  Table  5.  4.  1. 
These  rates,  together  with  the  flows  of  programs  and  results  shown  in 
Fig.  5.4.  1,  define  the  statistical  load  sharing  policy.  Note  that  due 
to  the  large  ratio  ofMp/Mr>  the  flows  of  programs  and  results  are  not 
symmetric.  In  this  case,  in  fact,  no  communication  link  carries 
both  programs  and  results.  The  solution  was  obtained  using  the 
algorithm  of  Section  5.3  and  is  guaranteed  to  be  within  1 % of  the 
optimum  load  sharing  policy.  This  solution  required  0.  3 sec  of 
computation  time  on  an  IBM  370-168  using  a FORTRAN  G1  complied 
program.  The  shortest  route  computation  times  were  kept  minimal 
by  the  use  of  the  decomposition  algorithm  of  Chapter  III. 

5 . 5 Additional  Remarks  on  the  Approximations  of  Section  5.2 

This  section  serves  two  principal  purposes.  The  first  is  to 
expand  on  the  rationale  underlying  the  approximations  of  Equations 
5.2.  1 and  5.2.2,  and  the  second  is  to  indicate  the  magnitude  of  the 
errors  due  to  these  approximations.  The  latter  is  accomplished 
through  several  mechanisms.  Namely,  the  error  in  5.2.1  is  derived, 
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the  nature  of  the  error  in  5.2.2  is  discussed,  and  the  changes  in  the 
assumptions  which  make  5.2.  1 and  5.2.2  exact  are  presented. 

In  Section  5.2,  Equation  5.2.  1 was  proposed  as  an  approxi- 
mation to  the  average  delay  (waiting  time  plus  service  time) 
experienced  by  a job  in  passing  through  the  ith  channel.  In  essence, 
one  postulated  two  independent  Poisson  input  streams  to  the  channels; 
(i)  the  program  stream  with  arrival  rate  f . and  service  rate  u C., 

pi  p i 

and  (ii)  the  result  stream  with  arrival  rate  f . and  service  rate  u C.. 

ri  'p  i 

Equivalently,  one  may  model  these  two  streams  as  one  Poisson 
stream. with  arrival  rate  (f  ^ + f^.)  and  a probability  density  function, 
pt(t  ).  of  the  service  time,  t,  for  a randomly  selected  arrival  given 
by: 


PtV 


-■P* 

f . + f . 

pi  ri 


f . 

M C.  exp(-M  C.t  ) + -s — 2L-  M C.  exp  (-n  C.t  ) 
pi  p i o fpi+  ri  r 1 r i o 

for  t ^0 
o 

otherwise 


One  can  see  that  the  assumptions  imply  that  the  channel  behaves  as  an 
M/G/l,  first-come-first-serve  (FCFS)  queueing  system. 

The  average  delay  through  this  system  is  derived  as  [KLEI76]: 


(fpi  + fri)E(t  > 

2(1  - (f  . + f JE(t)) 


P1 


rr 


E(t) 


or  equivalently 
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f . f . 
-Ei  + -£i 


f f ■> 

P'ri  M 

V "r 


T~ 


(f  . + f .)  (c.  - (-22. 

pi  ri  V i \n 


T 


ri  N\ 


(5.  5.  1 ) 


Upon  a comparison  of  5„  5.  1 to  5.2.1,  it  is  clear  that  5.2.1  is  a good 

approximation  to  5.  5. 1 when  the  third  term  in  the  numerator  of  5.  5.  1 

is  small  compared  to  the  sum  of  the  first  two  terms.  Although  that 

is  true  in  some  cases,  as  when  u =»  u , there  are  cases  in  which  the 

'p  r 

third  term  is  much  larger  than  the  sum  of  the  first  two  terms.  For 

example,  take  C.  = 1,  u = 0.  1,  f . = 0.04,  a =10,  and  f . = 4;  in 
r i r ri  P pi 

9 

this  case,  the  third  term  is  (0.  16)(9.  9)  - 15.68  and  the  sum  of  the 

first  two  terms  is  0.  4 + 0.  4 = 0.  8. 

Now  concentrating  on  Equation  5.2.2,  it  can  be  understood  to  be 
an  approximation  for  two  reasons.  The  first  contributor  to  the 
inexactness  of  5.  2 . 2 is  the  approximation  of  5.  5.  1 by  5.  2 . 1 , and  the 
second  is  the  imprecise  application  of  the  theory  of  networks  of  queues. 
From  an  equivalent  point  of  view',  the  ith  channel  is  assumed  to  behave 
as  an  M/M/1  queueing  system  with  service  rate  C.(f  , + f .)/( f .la  - 

l pi  ri  M pi  Pp 

fri/Mr)*  This  single  assumption  leads  to  the  approximation  designated 
5.2.  1,  and  it  also  permits  the  application  of  the  theory  of  networks  of 
M/M/1  queues  in  order  for  5.2.2.  to  result. 

Strictly  speaking,  the  form  of  5.2.2  is  founded  on  the  theory  of 
networks  of  queues  which  allows  treatment  of  each  queue  in  the  system 
as  if  it  were  independent  of  all  other  queues  in  the  system,  but  only  in 
the  case  that  each  queueing  center  has  a special  characteristic.  In 
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particular,  if  the  input  stream  to  any  given  queueing  center  is  a 
Poisson  process,  then  the  departure  process  must  also  be  Poisson 
[KOBA77,  MUNT73].  However,  a channel  behaving  as  an  M/G/l 
FCFS  queueing  system  is  known  not  to  have  this  property  [KLEI76]  , 
causing  difficulties  in  the  justification  of  5,2.2.  Therefore,  not  onlv 
are  certain  terms  in  5.2.2  inexact,  but  the  entire  form  of  5.2.2  is  an 
approximation,  so  that  there  would  seem  to  be  no  simple  manner  in 
which  to  evaluate  the  associated  error.  Of  course,  the  problems  of 
approximations  vanish  when  the  solution  has  a certain  property^,  i.  e. 
an  absence  of  channels  which  are  used  by  both  programs  and  results. 

On  the  other  hand,  one  could  have  altered  the  manner  in  which 
the  channel  operates,  and  thereby  avoided  all  of  these  issues.  In  such 
a case,  5.2.1  and  5.2.2  would  be  exact  results.  The  required 
change  is  to  operate  the  channels  in  a preemptive-resume  last-come- 
first-serve  (LCFS)  mode  or  in  a round  robin  (RR)  capacity  sharing 
mode  [BASK75].  In  actual  implementation,  both  techniques  require 
more  complex  nodal  processing.  On  the  other  hand,  there  are  some 
beneficial  effects.  By  giving  short  messages  a higher  priority,  one 
has  attained  generally  lower  expected  delays  (evidenced  by  the 
comparison  of  5.5.1  with  5.2.1).  Also,  the  analysis  and  predictions 
of  delays  are  exact  subject  onlyr  to  the  conditions  and  assumptions  of 
Section  5.  1 . 


Computer 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


X (jobs/sec) 

Processing  rate  of 
jobs  at  each  computer 
with  load  sharing 
(jobs  / sec ) 

. 00400 

. 01671 

. 00339 

. 01532 

. 00942 

1 

.01718 

.02428 

. 01814 

. 00594 

. 01847 

. 03490 

.01975 

. 04383 

.01949  1 

I 

. 02103 

. 01902 

. 00267 

. 01803 

.03111 

. 01847 

- 1 

Table  5.  4.  1 Computer  input  rates  for  example  shown  in 
Fig,  9.  The  expected  job  time  T for  this 
example  is  infinite  before  load  sharing  and 
230.2  sec/job  after  load  sharing. 


CHAPTER  VI 


CONCLUSIONS 

This  thesis  has  been  concerned  with  the  general  routing  problem 

in  computer-communication  networks.  One  of  the  best  known 

algorithms  which  solves  the  static  minimum  average  delay  problem  is 

the  EF  algorithm.  That  algorithm  has  been  modified  to  make  the 

routing  information  implicit  in  the  notation,  rather  than  requiring  the 

maintenance  of  an  explicit  routing  matrix  for  each  column  involved  in 

the  solution.  The  two  modifications  which  were  presented,  namely 

the  TR  and  CH  algorithms,  were  shown  to  be  more  compact  in  terms 

of  memory,  while  retaining  the  rapid  convergence  properties  of  the  EF 

algorithm.  Either  algorithm  would  seem  to  be  preferable  to  the  EF 

algorithm  for  most  networks;  generally  the  TR  algorithm  would  seem 

* 

to  be  preferable  to  the  CH  algorithm  for  sparse  networks,  while  the 
opposite  is  true  for  densely  connected  networks. 

Chapter  III  describes  two  new  shortest  route  algorithms,  the 
NXN  and  IHU  algorithms,  which  are  well  suited  to  the  subproblems  of 
the  various  routing  algorithms  of  this  thesis.  It  is  characteristic  of 
those  routing  algorithms  that  (i)  the  shortest  route  calculations  are 
done  many  times  for  the  same  network  topology,  and  (ii ) the  storage 
i of  the  entire  distance  matrix  is  not  burdensome  compared  to  the  other 

memory  requirements.  Thus  a modest  investment  in  the  decomposition 
of  the  network,  and  the  allocation  of  a small  amount  of  memory  to 
record  the  information,  can  result  in  substantial  reductions  in  shortest 


-119- 


. 1 •■>0- 


I route  computation  time.  The  N.\N  algor ithm  is  most  effective  for 

networks  which  are  planar  and  sparsely  connected,  and  the  IHC 
algorithm  is  most  effective  for  those  networks  of  a more  particular 
structure  described  in  Chapter  III. 

In  Chapter  IV,  a linear  routing  problem  was  defined  which  mini- 
mizes the  maximal  utilization  factor  of  any  channel  in  the  network. 

Such  a routing  policy  would  seem  desirable  in  a computer-communica- 
tion network,  especially  when  considering  the  finite  buffer  size  at  each 
node  of  the  network.  This  follows  from  the  observation  that  mini- 
mizing the  maximum  channel  utilization  corresponds,  approximately, 
to  minimizing  the  maximum  probability  that  a packet  is  lost  at  any- 
given  node  due  to  lack  of  buffer  facilities.  The  minimax  routing 
problem  was  solved  initially  as  a linear  program  via  decomposition 

i 

techniques  analogous  to  those  of  the  nonlinear  algorithms  of 
Chapter  II.  It  was  generally-  found  that  the  LCH  algorithm  was  the 
most  effective  on  those  problems  which  are  difficult  to  solve,  and 
hence  on  those  problems  where  efficiency  is  the  most  important. 

These  linear  algorithms  were  then  compared  to  a nonlinear 
algorithmic  approach  based  on  the  algorithpis  of  Chapter  II.  If  one  is 
content  with  the  inaccuracy  of  the  nonlinear  algorithms  (generally  less 
than  2%  of  the  objective)  at  the  ''knee"  of  the  convergence  curves,  then 

I 

they  are  preferable  to  the  linear  algorithms  as  the  complexity  of  the 
problem  increases.  In  particular,  the  nonlinear  algorithms  automat- 

I 

ically  take  partial  account  of  higher  order  criteria  for  the  solution; 
and  they  more  easily  solve  a perturbation  of  a problem,  using  the 
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I previous  estimate  of  the  solution  as  a starting  point,  than  the  linear 

programs  do. 

Finally,  in  Chapter  V,  an  example  of  the  application  of  the  ideas 
of  the  thesis  to  computer  load  sharing  was  developed.  Although  the 
problem  at  first  glance  may  not  seem  to  be  a multicommodity  flow 
problem  in  the  sense  of  Chapter  I,  it  can  be  modeled  as  such.  The 
example  serves  to  indicate  how,  for  a particular  problem,  the  special 
structure  of  that  problem  may  be  utilized  to  increase  the  efficiency  of 
the  algorithms  presented  in  this  thesis. 

Suggestions  for  future  research.  This  final  chapter  will 
conclude  with  selected  suggestions  for  future  research.  -First,  in 
Chapter  II,  the  range  of  algorithms  under  consideration  was  rather 
narrowly  constrained.  For  example  the  question  of  whether  to 
centralize  or  distribute  the  algorithms  for  rouring  was  not  asked,  and 
remains  unanswered.  To  define  the  "robustness"  of  an  algorithm  in 
such  a way  as  to  permit  a meaningful  comparison  of  distributed  and 

I 

centralized  algorithms  for  a fixed  allocation  of  the  resources  of  the 
network,  would  have  a great  impact  on  the  computer- communication 

network  area. 

Second,  there  may  be  ways  to  improve  the  convergence  of  the 
algorithms  of  Chapter  II  when  they  are  in  the  neighborhood  of  the 
solution.  Finding  downhill  directions  for  subgroups  of  commodities, 

1 rather  than  for  the  entire  network,  may  be  one  direction  of  research. 

Insofar  as  how  commodities  should  be  grouped,  the  only  limitations 
seem  to  be  those  imposed  by  one's  perception  of  the  cause  of  the 
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decrease  in  experimental  convergence  rate  of  the  various  algorithms. 

Most  extensions  of  work  in  Chapter  II  would  have  direct  applica- 
tion of  Chapter  IV.  But  Chapter  IV  has  some  problems  which  are 
unique.  In  particular,  improved  procedures  for  estimating  the 
bottleneck,  in  order  to  improve  the  lower  bound  of  Section  4.4,  would 
be  of  interest.  There  may  exist  characteristics  of  the  bottleneck 
channels  which  could  be  exploited  in  order  to  narrow  the  choice  of 
channels  to  be  included  in  the  estimate. 

Appendix  B mentions  "stability"  of  a network  without  attempting 
a definition  of  the  term.  To  define  stability  and  to  pursue  its 
implications  could  have  great  impact  on  the  problem  of  dvnamic  net- 
work control.  Depending  on  the  direction  which  such  a study  might 
take,  some  of  the  algorithms  of  this  thesis  might  have  great  utility  in 
finding  the  routing  policies  which  maximize  the  stability  of  a network. 


APPENDIX  A 
SELECTED  PROOFS 

This  appendix  contains  the  proofs  of  selected  remarks  within  the 
body  of  the  report  where  continuity  of  thought  prevented  a more 
rigorous  discussion. 

u 

Theorem  2.  1 . A:  The  absolute  delay  error,  T(f'  ) - T(i  ),  at  the  com- 

pletion of  a restricted  master  problem  is  upper  bounded  by 
- Rel  Cost  where  f is  the  current  estimate  of  f , fff  is  the 

solution  to  2.  1.1,  and  is  the  extremal  flow  generated  by  the 

associated  subproblem. 

Proof:  Since  T(f)  is  a convex  function,  then  it  lies  above  the  best 
linear  approximation  of  T(f),  i.  e. 

T(f)  ^ T(f*)  + <VT(f*),  f-f*> 

If  f<F,  then 

T(f~)  - T(f)  £ - <VT(f*),f-f*>  £ - Rel  Cost  (<p(;j+1)) 

U 

by  the  definition  of  Rel  Cost  In  particular,  let  f = f , which 

proves  the  theorem. 

Theorem  2. 2.  A;  Any  flow  f which  can  be  expressed  as  a convex  com- 
bination of  trees,  can  be  expressed  as  a convex  combination  of  at  most 
NA  + N of  them. 

Proof:  If  f can  be  expressed  as  a convex  combination  of  trees,  then 
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there  exists  tX-  } such  that  f = L X,i„  <A,k>,  S X , = 1,  and 
ik  k = 1 i = 1 lk  (i)  i = i Aik 

Xlk  2 0,  where  the  notation  is  in  agreement  with  that  of  Section  2.2. 
Using  linear  programming  terminology,  [\lki  is  a feasible  solution  to 
a system  of  NA  + X linear  equations  which  are  shown  in  Equation  4.  1. 3; 
therefore,  the  system  Iiu.j  u.  basic  feasible  solution,  (x1.,  },  where  there 

1 K 

are  at  most  XA  + N nonzero  elements  [DANT63],  (Thus,  if  the 
number  of  trees  in  the  basis  of  the  TR  algorithm  exceeds  NA  + N,  then 
the  basis  may  be  reduced  to  NA  + N elements  by  pivot  operations  which 
are  fundamental  in  linear  programming.  ) 


Theorem  2,3.  A:  The  direction,  Af,  as  defined  in  section  2.3,  is  a 

feasible  downhill  direction. 


* * N Jk  k 

Proof:  Recall  the  form  of  f , i.e.  f = .L.  X.,  ft  , such  that 

jk  k = 1 i = l ik  (i) 

xik  Any  flow,  f,  which  can  be  expressed  in  the  same 

form  (as  a convex  combination  of  the  trees  in  the  basis),  is  said  to  be 

feasible  with  respect  to  the  current  basis.  Af  is  a feasible  direction 

if  there  exists  & > 0 such  that  f + fiAf  is  feasible  with  respect  to  the 

current  basis.  Assume,  momentarily,  that  everv  >0;  then  bv 
* N tk  k 

construction,  f + fiAf  = .2^  ( X-k  + 6AX^)<P^j.  From  the  definition 

Axik  = -8T(f")/axik+  U/y  J^OTfD/ax^,  the  condition  .£  (Xik  + 

6AX,„)  = 1 is  maintained  for  any  5 . Also,  if  0 s ^ s min 

Ax,  < 0 lk  lk 
Aik 

then  (Xik  + 6Axik)  s 0.  Thus,  Af  is  a feasible  direction.  Notice  that 
the  restriction  X-k  > 0 may  be  relaxed  to  X^k  2 0 if  AX-k  > 0.  Thus  a 
newly  added  from  a subproblem  poses  no  difficulty  even  though 
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\ik  = 0,  since  Axik  ■>  0. 

Af  is  a downhill  direction  if  (f*  + 0Af>  fi=  Q < 0.  Furthermore. 


a?(f  + /3Af)^  = 0 = <VT(f"),Af)  = S=i  .2  AXjk 


= E E + J_  g &T(f*)N  6T(f*) 

k=li--ls  6\ik  tkm  = l 6xmk-  5Xik 


N 
= E 
k= 


g / k . m£if  + 1 / £ p£if } * 0 

= 1 *-  i=l  ' ^ ^-llr  ' tlc  s i=  1 *Xil,  / / 


where  the  last  inequality  follows  from  the  Schwartz  inequality. 

Equality  occurs  only  if  tfT(f*)/ax1>k  = 0T(f^)/£>X2#  k = •••  = £T(f*)/6Xt  k 

for  every  kc£l,  2.  • * • , N) , which,  in  turn,  implies  f*  is  the  solution  to 
the  current  master  problem  ot  the  TR  algorithm  in  which  case  there 
can  be  no  feasible  downhill  direction.  Q.  E.  D. 

Theorem  2.  4.  A;  Any  flow  f which  can  be  expressed  as  a convex  com- 
bination of  chains,  can  be  expressed  as  a convex  combination  of  at 
most  NA+N(N- 1 ) of  them. 


Proof;  If  f can  be  expressed  as  a convex  combination  of  chains,  then 

there  exists  {x.  } such  that  f = 8 S ^ \ <pmk  C£ k \ -1 

imk  m=l  k=l  i = l umk  ti)  Ximk  ** 

k ^ m 

and  ximk  ^ 0 where  the  notation  is  in  agreement  with  that  of  Section 
2.4.  The  proof  follows  the  same  lines  as  Theorem  2. 2.  A after  one 
notes  that  the  set,  tX]mk}.  is  a feasible  solution  to  a system  of  NA+(N-1)N 
linear  equations  analogous  to  those  shown  in  Equation  4.  1.  3. 
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Theorem  4.  4.  A:  The  lower  bound,  S,a,  f)  as  defined  in  Section 

4.  4,  can  converge  to  $ . 

Proof:  Whether  or  not  a,  f)  converges  to  depends  on  the 

methodology  employed  in  the  selection  S and  a of  ^(S ,a , f)  and  e of 
AL*(a,  e,y).  In  this  proof,  it  is  shown  that  there  exists  at  least  one 
set  of  choices  which  implies  the  convergence  of (S,u,  f)  to  $ . 

Assume  that  the  lower  bound  is  computed  upon  the  completion  of 
step  2 of  the  AP:(a,  e,  y)  algorithm  on  each  iteration,  and  that 
(1)  a • aK 

(ii)  f . fK 

(iii)  S is  chosen  as  the  set  of  all  the  channels  to  the  network. 

(iv)  The  parameter  e is  chosen  such  that  0 < c < 1. 

IS 

Under  these  circumstances,  the  lower  bound  is  denoted  , and 


2;  (j 

PL  K 


N A a , . 
K 


NA 

(1  - c)  L 


i = 1 1 - p./o 

i K 


* o 


NA  a,, 
h 


K 


(1  - €) 


1 - $Y(f^)/aK 


* a 


K 


rrh)  taK  - <VfK» 


Now  it  is  known  from  the  proof  of  the  convergence  of  the  AL*(a,  C,y) 
algorithm  (Equation  4.  4.  1)  that; 


where  r = a /(NA(1  + e)).  Since  aK  - ^ (f*)  s aR  - then 

a ^1  ' (I  ‘ r)K  1 ^(fu)  ) 

The  last  inequality  and  the  fact  ^ ^ imply  that  0*  converges  at 

least  geometrically  to  $ . Q.  E.  D. 


APPENDIX  B 

PERTURBATIONS  AND  THE  ROUTING  ALGORITHMS 

This  appendix  demonstrates  the  capability  of  routing  algorithms 
to  resolve  the  minimum  average  delay  problem  for  a perturbation  of  a 
demand  matrix  for  which  the  problem  has  already  been  solved.  As 
noted  in  Chapter  II,  such  capability  is  inherent  in  the  EF,  TR,  and  CII 
algorithms  although  the  FW  algorithm  does  not  have  this  ability. 

There  is  no  reason  to  believe  that  the  EF,  TR  and  CH  algorithms 
should  behave  very  differently  from  une  another  on  a perturbation  of  a 
previously  solved  problem;  hence,  for  convenience,  the  methodology 
for  handling  perturbations  (as  well  as  the  example  of  its  use)  is 
presented  only  for  the  TR  algorithm. 

As  a by-product  of  the  investigation  of  perturbations  of  the 
requirement  matrix,  one  can  obtain  a qualitative  feel  for  the  proper 
time  intervals  between  updates  of  the  routing  policy  for  a quasistatic 
implementation.  This,  however,  would  need  to  be  conditioned  on  the 
probabilistic  nature  of  the  elements  of  the  requirement  matrix,  and 
only  one  model  is  considered  here.  The  reasons  that  the  particular 
model  was  chosen  may  be  more  clear  after  reading  Appendix  C. 

A model  of  the  requirement  matrix.  As  a preliminary,  the 
model  of  the  perturbations  of  the  demand  matrix  is  presented.  The 
requirement  matrix  will  be  defined  at  discrete  times,  tc{0,  1 , 2,  • • • } , 
and  will  be  denoted  R(t)  = fr_(t)]  . The  following  specify  the  demands. 
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i Q n 
i — J 


(i)  r-j(tj)  and  , [n^(t9)  are  independent  random  variables 

• providing  1 * m or  j / k for  every  t and  t9. 

1 2 

(ii)  If  i = j,  r.j(t)  = 0.  Otherwise,  r^.(t)  has  the  density 
function  of  the  absolute  value  of  a zero  mean 
Gaussian  random  variable  with  a variance  of  3. 


(iii)  If  i / j,  then  r (t)  given  r (t  - 1)  has  the  density 

J A J 

function  of  |r.  (t  - 1)  + a x (t)  1/(1  i-  a2)% , where  x (t) 

XJ  1 t 1J 

is  a zero  mean  Gaussian  random  variable  with  a 
variance  of  3 which  is  independent  of  all  other  random 
variables  in  the  network  (including  x.^lT)  if  r / t,  and 
where  at  > 0.  It  is  not  hard  to  verify  that  (ii)  and 
(iii)  are  not  in  conflict. 


The  problem.  The  problem  which  is  to  be  solved  is  described 
below.  The  network  is  the  seven  node  symmetric  network  of  Fig. 
4.3.3.  Every  channel  is  assumed  to  have  capacity  5.5.  The  require- 
ment matrix,  R(t)  = [r  (t)],  at  t = 0 is  the  one  which  would  be 
specified  by  the  rules  of  Appendix  C for  a static  problem.  The 
minimum  delay  problem 


min 

f€F 


f. 

i 


5.  5 - f. 

l 


is  solved  such  that  the  upper  bound  on  the  error  is  less  than  2%.  Then 
R(l)  is  generated  according  to  specification  (iii)  above.  The  current 

N \ . 

estimate  of  the  solution,  i.e.  f';:(0)  = L l \ 

k=  1 j = l Jk 


L 


is  modified  in 


* 


the  following  way:  the  X remain  fixed,  and  the  o'  are  changed  - 

lK 

R(l)  were  the  requirement  matrix  when  they  were  generated.  Th.~ 
similar  (but  not  equivalent)  to  keeping  the  routing  policy  fixed  w.-ule 
the  requirement  matrix  has  changed. 

Clearlv.  th>-  new  ••••-mlting  flow,  f/ll,  is  a member  of  ♦!;  r.  " 
set  F;  however,  the  flow  may  not  be  feasible,  in  which  case 
appropriate  measures  are  taken  (as  indicated  in  Section  4.  5).  Since 
the  mapping  of  the  requirement  matrix  into  a flow  for  a fixed  set  of 
parameters,  },  is  continuous,  if  R(l)  is  ’’close"  to  R(0),  then  f(  1 ) 
is  "close"  to  f (0),  and  f(l)  is  unlikely  to  be  unfeasible.  The  new 
minimum  delay  problem  is  solved  to  a 2 ,r<  tolerance,  and  the  procedure 
is  iterated  for  R(2),  R(3),  R(4),  and  R(5>. 

Table  B.  1 summarizes  the  convergence  of  the  TR  to  the  succes- 
sive 2%  tolerance  estimates  of  the  solutions  to  the  given  minimum 
delay  problem.  The  following  points  are  noted.  The  load  on  the  net- 
work is  fairly  heavy  with  some  channels  in  each  solution  having  a 
utilization  factor  greater  than  0.  7.  Resolving  the  problem  for  a per- 
turbation is  generally  fast  compared  to  solving  the  problem  from  an 
arbitrary  starting  point  depending  on  the  size  of  the  perturbation. 

For  this  problem,  the  degradation  of  the  expected  delay  from  the 
optimum  value  would  generally  be  less  than  4 or  5 °o  for  an  optimum 
routing  policy  computed  1 second  beforehand  if  £ 0.  10  and,  at  least 
to  the  author,  this  was  surprisingly  low.  It  is  also  noted  that  the 
random  variables  generated  for  the  problem  were  not  exactly  Gaussian, 
and  one  effect  of  this  was  a drift  towards  a more  heavilv  loaded  network 
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on  each  successive  minimum  delay  problem. 

A parenthetical  remark.  As  one  looks  at  the  above  problem, 
naturally  the  question  of  the  sensitivity  of  a routing  policy  to  changes 
in  requirements  or  inaccurate  estimates  arises.  If  f solves  B.  1,  it 
is  well  known  that  the  routing  policy  which  corresponds  to  f is  not 
unique.  If  one  says  that  he  wants  the  "best"  of  the  many  routing 
policies,  then  a dilemma  arises  as  to  how  the  best  should  be  defined. 

On  the  one  hand,  the  routing  policy  with  the  greatest  number  of  "pure" 
routing  strategies  has  the  advantage  of  reducing  out-of-order  message 
deliveries.  That  is,  if  the  packets  belonging  to  a given  message 
traverse  different  routes,  then  their  arrival  at  the  destination  is 
generally  out  of  sequence,  which  complicates  the  reconstruction  of  the 
message.  However,  it  is  intuitively  clear  that  the  routing  policy 
which  is  the  most  split  is  also  the  least  sensitive  to  arbitrary  perturba- 
tions of  the  requirement  matrix.  This  results  in  increased  intervals 
between  necessary  routing  updates,  which  may  be  interpreted  as 
increasing  the  stability  of  the  network. 

The  intuitive  clarity  of  this  is  enhanced  by  the  following  argument. 

Let  each  r..  be  an  independent  random  variable  with  E(r..)  = r and 
i]  U 

2 

Var(r.j)  = a . On  the  ith  channel,  perhaps  one  is  given  the  choice  of 

routing  commodity  (1,2)  alone,  or  the  sum  of  the  fractions,  1 '(X-l), 

of  each  commodity  (1,  k)  for  ke  {2,  3,  • • • , n}  . The  expected  value  of 

flow  on  the  ith  channel  is  the  same  in  each  case,  but  the  variance  of  the 
2 2 

flow  is  a for  the  former  and  CT  /(N-l)  for  the  latter.  It  follows  that  if 
the  same  types  of  choices  go  on  throughout  the  entire  network,  and  one 
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alvvavs  chooses  to  split  routes  as  much  as  possible,  then  one  has  in 
some  sense  minimized  the  variation  of  the  flow  in  the  network  from 
that  which  is  expected.  This,  in  turn,  increases  the  stability  of  the 
network. 


Value  of  the 

Value  of  the 

Computation 
time  (seconds! 
required  for  the 

expected  delay 

expected  delay 

new  estimate  of 

using  the  set. 

after  minimiza- 

the solution  of 

U-nJ.of  the 

tion  to  a maxi- 

the  minimum 

l 

estimate  ui  the 

m cl  1 2 o toler- 

delay  problem  j 

solution  to  the 

ance  in  the 

using  the  previ- 

previous 

error. 

ous  set,  { Xjk-  . 

problem. 

as  a starting 
point  for  the 
minimization. 

0. 63546 
0. 63949 
0. 64533 
0. 64705 
0. 66373 
0. 73243 


0.  55 
0.  06 
0.  09 
0.  07 
0.  09 


Table  B.  1 A summary  of  the  results  of  the  perturbation  problem 
on  the  seven  node  symmetric  network  using  the  TR 
algorithm. 


I 


APPEMJIX  C 


THE  REQUIREMENT  MATRICES  I -'OR  THE  SAMPLE  PROBLEMS 


This  appendix  specifies  the  requirement  matrices  for  the 
problems  of  Chapter  IV.  Some  choices  made  in  this  appendix  are 
arbitrary.  There  is  an  attempt  here  to  explain  how  they  were 
and  to  justify  them. 

The  first  consideration  was  to  avoid  extraordinary  degeneracies 
that  might  arise  in  solutions  to  4.  1.  1.  Such  degeneracies  would  arise 
if  one  made  a typical  assumption  that  all  off-diagonal  elements  or  the 
requirement  matrix  are  equal.  Furthermore,  the  assumption  is 
counter-intuitive  when  compared  to  the  assumption  that  all  such 
elements  merely  have  the  same  statistics.  The  latter  assumpt.cn 
was  employed. 

In  choosing  the  statistics  of  the  elements,  convenience  was  a 
primary  objective.  A computer  program  devised  by  Rader  {RADE6P) 
which  generates  uncorrelated,  nearly  Gaussian  pseudo- random- 
variables  met  this  criterion.  This  allowed  regeneration  of  the 
matrices  for  different  algorithms  that  solved  the  same  problems 
without  requiring  maintenance  of  a large  table  of  random  varia:  U s. 
Moreover,  a simple  method  exists  for  perturbing  the  matrices  ir.  a 
quite  arbitrary  fashion  while  retaining  the  form  of  the  density  functions 
of  the  demands  (see  Appendix  R). 

In  keeping  with  the  above,  each  term  r of  the  demand  matrix 
was  assumed  to  be  a random  variable  with  the  following  properties: 
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(1)  r.  . and  r , no  independent  random  variables  in  the 

ij  mk 


case  i t m or  j r k. 


Oil 


r . . = 0 if  i 
’-J 


Otherwise,  the  probability  densitv 

of  r.  . is  given  bv: 
hi  S 

’ r> 

l (2/ 371)^  exp  (-r“/6)  for  r -i  0 
1 r o o 


Pr.^o1  - | . . 

i]  ( 0 otherwise 


That  is,  r . . x. . if  i ^ i,  where  x.  . is  a zero 
ij  ij  ij 


mean  Gaussian  random  variable  with  a variance  of  3. 
The  results  of  Chapter  IV  are  independent  of  the  value 
of  the  variance  except  as  a scaling  factor. 


Clearly,  these  properties  of  the  elements  of  the  requirement 
matrix  are  not  suited  to  every  computer- communication  network. 
Among  other  possibilities,  a probabilistic  model  of  the  elements  for  a 
particular  network  might  include:  (i)  dependency  among  the  rancom 

variables,  (ii)  a nonzero  probability  that  r = 0 of  i * 3,  (iii)  an  unper 
bound  on  the  elements,  and  (iv)  different  means,  variances,  forms  of 
densities,  etc.,  which  depend  on  i and  j.  Nevertheless,  the  author  is 
not  of  the  opinion  that  the  algorithms  of  this  thesis  are  highly  dependent 
on  the  model  of  the  requirement  matrix. 

Table  C.  1 displays  a set  of  nearly  independent  samples,  { v } 

m k 


(where  y , is  the  kth  element  of  the  mth  row  of  the  table),  of  a random 
J mk 


variable  y which  is  zero-mean,  approximately  Gaussian  with  a variance 
of  3 as  generated  by  Rader's  algorithm.  The  requirement  matrices 
of  the  sample  problems  in  Chapter  IV  are  specified  by: 


J 


m = |(j  • (i  - 1 )\)/8f 


and 


k = .1  + (1  - 1)N  - 8(m  - 1) 


where  [x]  is  the  smallest  integer  which  is  not  less  than  x.  The 
definitions  of  m and  k are  simply  specifications  of  the  mapping  of  the 
table  of  random  variables  into  the  demand  matrix. 

For  example,  for  the  seven  node  symmetric  network,  the  value 
r^_  is  desired.  Thus, 


m 


f(5  + (3  - 1 )?)/ 8 ] 
19-8(3-1)  = 


[2.875! 


which  vields  r„^ 
3 a 


l y33  ' 


0.  11183. 
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Table  C.  1 These  values  specify  the  elements  of  the  requirement 
matrix  for  the  sample  problems  of  Chapter  IV 
according  to  the  rules  of  Appendix  C. 
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A I* PEND IX  D 


DETERMINATION  OF  R(  >1  I ING  FOR  HIGHER  ORDER  CRITERIA 

This  appendix  presents  one  example  of  a linear  routing  problem 
solved  for  higher  order  criteria.  The  results  which  are  presented 
are  typical  of  those  problems  for  which  higher  order  criteria  were 
considered.  The  example  is  of  the  ten  node  problem  considered  in 
Section  4.3  and  redrawn  in  Fig.  D.  1.  In  that  figure,  along  each 
arrowed  arc  is  a pair  (i,  /3 ) where  i is  the  first  criterion  for  which  the 
respective  channel  was  associated  with  a negative  dual  variable,  and 
ft  is  the  value  of  pL  The  problem  was  solved  through  the  fourth 
criterion. 

Table  D.  2 presents  the  computation  times  for  the  various 
criteria  with  the  LTR  algorithm.  Each  iteration  of  the  linear 
programming  results  which  is  listed  represents  one  of  two  types  of 
events: 

(i)  If  the  iteration  is  within  the  same  criterion  level 
as  the  previous  iteration,  columns  have  been 
added  and  the  linear  program  has  been  remini- 
mized.  For  example,  on  the  5th  iteration  ten 
columns  (which  in  this  case  are  trees)  were  added 
to  the  program,  and  the  best  value  of  the  objective, 
subject  to  the  31  columns  presently  available,  was 
determined  as  10.48. 

(ii)  If  the  iteration  is  on  a different  criterion  level 
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Fig.  D.  1 A ten  node  network  with  thirty  directed  arcs. 


Criterion,  i | Iteration 


Value  of  the 

OujOCtl  i',  H- 


Cumulative 

computation 


Table  D.  1 Computation  time  for  the  LTR  algorithm  to  solve  the 
linear  routing  problem  associated  with  Fig.  D.  1 
through  the  fourth  criterion. 
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from  the  previous  iteration,  the  program  was 
changed  to  reflect  the  new  criterion,  and  the  new 
objective  minimized  with  the  columns  existing 
from  the  previous  iteration. 

It  is  noted  that  the  second  type  of  iteration  is  particularly  time 
consuming.  This  reflects  the  large  number  of  basis  changes  required 
to  accomplish  the  minimization.  Moreover,  as  the  number  of  the 
criterion  increases,  the  number  of  active  constraints  (which  is  directly 
related  to  the  size  of  the  basis  inverse)  generally  increases  as  well. 

This  is  because  every  channel  which  had  its  utilization  factor  fixed  by 
a previous  criterion  corresponds  to  a necessarily  tight  constraint  for 
the  current  program.  The  consequence  is  that  each  successive 
criterion  generally  requires  more  computation  time  than  its  predecessor. 

Also  in  Fig.  D.  1,  there  corresponds  to  many  of  the  channels  a 
pair  [x,y]  where: 

x = the  value  of  the  flow  in  the  channel  at  the  conclusion 
of  the  6th  iteration  of  the  LTR  algorithm, 
y = the  value  of  the  flow  in  the  channel  at  the  conclusion 
of  the  TR"  algorithm  (i.e.  after  1 second  of  compu- 
tation time)  which  produced  the  appropriate  curve 
in  Fig.  4.  5.  1 . 

The  pairs  are  provided  for  those  channels  for  which  either  x or  y is 
greater  than  9.  Observe  the  following:  for  each  set  {(i,p.),  [x,  y] } 

along  each  channel  in  Fig.  D.  1,  |y-j^|  s |x-p-  |.  This  verifies  the 
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suggestion  at  the  conclusion  of  Section  4.  5 that  the  AL  (a,  e,  y) 
algorithm  tends  to  produce  a solution  that  is  better  from  the  point  of 
view  of  higher  order  criteria  than  do  linear  programs  which  are 
solved  only  for  the  first  criterion. 


APPENDIX  E 


FURTHER  COMMENTS  ON  THE  LOWER  BOUND  FOR  THE 
LINEAR  ROT  TING  PROBLEM 


It  is  noted  in  the  discussion  of  Fig.  4.5.2  that  the  respective 
lower  bounds  for  the  TR  ',  EF*,  and  CII  algorithms  are  not  as  satis- 
factory as  those  of  the  other  examples.  This  is  true  even  though  the 
lower  bound  can  become  exact.  The  underlying  reasons  for  this 
phenomenon  are  discussed  in  Sections  4.4  and  4.  5.  This  appendix 
elaborates  on  that  discussion  and  demonstrates  the  possibility  of 
algorithms  that  do  somewhat  better. 

The  TR':  algorithm  was  run  for  1. 08  seconds  on  the  linear 
routing  problem  associated  with  Fig.  4.3.2,  at  which  point  ^ was 
91.73,  and  the  values  of  the  eleven  most  utilized  channels  were  as 


follows; 


Channel 

Utilization 

(20,  25) 

91. 73 

(24,  4) 

91.  58 

(14,  15) 

87.  43 

(25,  20 

86.  91 

(15,  16) 

86.  73 

(20,  22) 

86,  15 

(22, 20) 

86.09 

(23,  13) 

85.  33 

(1,2) 

83.  80 

(4,24) 

83.  05 

(15,  14) 

82.61 

To  look  at  this  list  and  conclude,  simply  from  the  values  of  utilization. 


-146- 


-147- 


that  the  top  three  channels  are-  the  bottleneck  would  be  absurd. 
However,  if  the  reader  refers  to  Fig.  E.  1 and  indicates  the  above 
channels  on  the  graph  in  the  order  of  appearance  in  the  table,  the 
following  observations  are  clear; 

(i)  The  first  cutset  occurs  after  the  third  channel  is 
considered.  One  would  therefore  say  that  this  is 
the  first  estimate  of  the  bottleneck  with  any 
potential  for  producing  a good  lower  bound  on  fi  . 

(ii)  The  addition  of  channel  (25,  20)  b}r  itself  will 

degrade  the  bound.  This  is  because  there  would 
also  be  a path  with  zero  estimated  bottleneck 
channels  between  nodes  25  and  20,  i.  e.  the  addition 
of  channel  (25,20)  could  not  form  or  contribute  to 
the  bottleneck. 

(iii)  The  addition  of  channel  (15,  16)  to  the  first  three 
(or  the  first  four)  channels  will  not  improve 
/S^(S,°°,f).  Two  channels  in  a simple  sequence  like 
(14,  15)  and  (15,  16)  are  not  beneficial  for  the  reasons 
discussed  in  Section  4.  4 

The  list  of  observations  can  go  on;  however,  it  seems  to  become 
rapidly  evident  that  the  three  most  utilized  channels  are  most  probab 
the  bottleneck,  and  the  calculation  for  the  lower  bound  based  on  any 
other  estimate  of  the  bottleneck  is  probably  in  vain. 
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To  put  tins  kind  of  subjective  a posteriori  intelligence  into  a 
computer  program  which  estimates  the  bottleneck  would  be  of  great 
benefit,  if  its  execution  required  little  computation  time  compared  to  a 
shortest  route  computation.  The  method  which  the  author  instructed 
the  TR*  algorithm  to  use  for  choosing  S,  for  the  previously  described 
run,  was  not  very  sophisticated.  It  was  as  follows;  (i)  choose  to  be 
all  channels  with  a utilization  factor  greater  than  (0.  9)^,  (ii)  calculate 
a lower  bound,  (iii)  reduce  the  cardinality  of  S by  a factor  of  0.  75 
(rounding  to  the  next  lower  integer  if  required)  by  rejecting  the  least 
utilized  channels,  and  (iv)  if  less  than  5 lower  bounds  have  been 
. calculated,  then  return  to  the  second  step.  The  calculated  lower 

I 
I 

bounds  were: 

' ' [Sj  BL(S,»,f) 

11  57.86 

8 44.23 

6 49.98 

4 66.84 

3 89.24 

The  last  of  these  values  is  in  fact 

J 

Although  this  appendix  probably  asks  more  questions  than  it 
answers,  it  is  meant  to  suggest  that  there  is  a great  deal  of  information 
on  the  identity  of  the  bottleneck  when  the  AL;:(a,  e,y)  algorithm  reaches 
the  knee  of  its  convergence.  Much  more  information  is  available  than 
is  used  to  produce  the  lower  bounds  in  Figs.  4.5.1,  4.5.2,  4.5.3,  and 
| 4.5.4.  The  problem,  of  course,  is  how  to  extract  that  information  in 

| 

an  efficient  manner. 
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